# Semi-analytical pricing of barrier options in the time-dependent
$\lambda$-SABR model.

Andrey Itkin 1 and Dmitry Muravey 2

1 Tandon School of Engineering, New York University, 1 Metro Tech Center, 10th
floor, Brooklyn NY 11201, USA
2 Moscow State University, Moscow, Russia

  * We extend the approach of Carr, Itkin and Muravey, 2021 for getting semi-analytical prices of barrier options for the time-dependent Heston model with time-dependent barriers by applying it to the so-called $\lambda$-SABR stochastic volatility model. In doing so we modify the general integral transform method (see Itkin, Lipton, Muravey, Generalized integral transforms in mathematical finance, World Scientific, 2021) and deliver solution of this problem in the form of Fourier-Bessel series. The weights of this series solve a linear mixed Volterra-Fredholm equation (LMVF) of the second kind also derived in the paper. Numerical examples illustrate speed and accuracy of our method which are comparable with those of the finite-difference approach at small maturities and outperform them at high maturities even by using a simplistic implementation of the RBF method for solving the LMVF.

††footnotetext: Keywords: barrier options, stochastic volatility, SABR model,
GIT method, semi-analytical solution, mixed linear Volterra-Fredholm equation,
radial basis functions††footnotetext: 2010 MSC: 91G40; 91G60; 91G80; 47G30;
35R09; 65L12;

## Introduction

The SABR model is one of the most popular models in mathematical finance.
Introduced in (Hagan et al., 2002) (and then an arbitrage-free version in
(Hagan et al., 2014)), it quickly became a standard tool among practitioners
especially in the interest rate derivative markets due to its ability to
capture skew and smile features observed in the interest rates implied
volatilities. For European options it has an asymptotic solution, (Hagan et
al., 2002) which allows closed form representation of the implied volatility.
Later a SABR extension for negative interest rates (the shifted SABR model)
was developed as well as advanced analytics was provided by using a heat
kernel expansion technique, see (Antonov et al., 2019) and references therein.
As mentioned in that book, since the SABR model allows an analytic
approximation for Black volatilities it is widely used for the swaption
volatility cube interpolation and extrapolation.

The SABR model can also be extended to its time-dependent version where all
model parameters are deterministic functions of time. This obviously
complicates the calibration procedure. An advanced calibration method of the
time-dependent SABR model based on so-called "effective parameters" was
developed in (Van der Stoep et al., 2015).

Another extension was proposed in (Hagan et al., 2020). The authors consider
the mean-reverting SABR model where the stochastic process for the
instantaneous volatility has a mean-reverting drift. They also assume that
coefficients of the model are functions of time. Then they develop asymptotic
methods to obtain an effective forward equation for the marginal density of
the forward price. This equation is not exact, but as accurate as the SABR
implied volatility formulas. The authors suggest solving this one-dimensional
partial differential equation numerically to obtain the density. European
option prices can then be found by integrating this density with the option
payoff.

As far as exotic options is concerned, for the standard SABR model with
constant coefficients in (Yang et al., 2017) a closed form expression is
obtained to approximate prices of various types of barrier option (Down-and-
Out/In, Up-and-Out/In). The authors derive an approximate formula for the
survival density and then represent the barrier option price as the one-
dimensional integral of its payoff function and the survival density. The
approximation error of the survival density is also analyzed.

Since even for the SABR model with constant coefficients the closed form
expression for the barrier option price is not available (without
approximations), various asymptotic methods have also been proposed assuming
that the model contains a small parameter, see e.g., (Barger and Lorig, 2017;
Kato et al., 2013) among others.

Alternatively, numerical methods have been used extensively to price options
under the SABR model especially when those options have high maturities and,
hence, the asymptotic solution of (Hagan et al., 2002) (for the European
option) becomes inaccurate. For instance, in (Thakoor et al., 2019) a
computational method based on a spectral discretization of the pricing
equation is developed for pricing options with discrete barriers under the
arbitrage-free SABR model. The high accuracy of the method is established by
comparison with special cases of the SABR model where analytical solutions are
available. But basically, in the literature various flavors of Monte Carlo and
finite-difference (FD) methods were used to price barrier options under the
SABR model with both time-dependent and constant coefficients.

As mentioned in (Hagan et al., 2020), the SABR model is effective at managing
volatility smiles, volatility as a function of the strike $K$ at a single
expiry date $T$. For each $T$, the SABR parameters $\alpha,\rho,\nu$ are
calibrated so that the model’s implied volatility curve
$\sigma_{I}\left(T,K\right)$ matches the market’s implied volatilities for
that $T$. However, the SABR model is less useful for managing volatility
surfaces, volatility as a function of the strike $K$ at multiple exercise
dates $T$, which are needed for handling many exotics. Managing volatility
surfaces requires a richer model, such as the dynamic SABR model where all the
model coefficients are functions of time.

Therefore, in this paper our main interest is pricing barrier options under
the time-dependent SABR model where the barrier level could also be a
deterministic function of time. We propose a slightly modified version of the
SABR model with the following changes: a) for the instantaneous stochastic
volatility we assume a lognormal vol-of-vol (as in the original SABR model)
but also add a mean-reverting drift; b) this drift has zero mean-reversion
level and is linear in $\sigma_{t}$. The form of the mean-reversion term for
$\sigma_{t}$ is chosen by a tractability argument. It is inspired by the
$\lambda$-SABR model of (Henry-Labordere, 2005). The assumption about the
linear drift can be relaxed, for instance a quadratic drift can be used
inspired by the so-called "3/2" stochastic volatility model, (Platen, 1997;
Carr and Sun, 2007)111In contrast to these papers all coefficients of our
model are arbitrary functions of time (in (Carr and Sun, 2007) the mean
reversion rate is an arbitrary function of time and the mean-reversion level,
vol-of-vol and correlation are constant). However, as this will be seen in
section 2 the linear drift allows reduction of 2D integrals to 1D ones, so the
corresponding LMVF equation becomes easier to solve.

With all these preliminaries in mind, we define the model by the following
stochastic differential equations (SDEs)222Note, that the SABR model has been
employed for Fixed Income, FX and Equities. Here without loss of generality we
formulate it for the Equity world, while, e.g. for FX this setting can be
found in (Van der Stoep et al., 2015).

$\displaystyle dF_{t}$ $\displaystyle=\sigma_{t}F_{t}^{\beta+1}dW^{(1)}_{t}$
(1) $\displaystyle d\sigma_{t}$
$\displaystyle=-\kappa(t)\sigma_{t}dt+\gamma(t)\sigma_{t}dW^{(2)}_{t},$
$\displaystyle d\langle W^{(1)}_{t},W^{(2)}_{t}\rangle$
$\displaystyle=\rho(t)dt,\quad[F,t]\in[0,\infty)\times[0,\infty),\quad
F_{0}=F,\quad\sigma_{0}=\sigma,$

where $t\geq 0$ is the time, $F_{t}$ is the stochastic forward price,
$\sigma_{t}$ is the stochastic volatility, $\kappa(t)$ is the rate of mean-
reversion, $\gamma(t)(t)$ is the volatility of volatility (vol-of-vol),
$\theta(t)$ is the mean-reversion level (the long-term run), $W^{(1)}$ and
$W^{(2)}$ are two standard correlated Brownian motions with the correlation
coefficient $\rho(t)\in[-1,1]$, and $\beta$ is the elasticity parameter such
that $-1<\beta<1,\ \beta\neq\\{0,-1\\}$333In case $\beta=0$ this model is the
Black-Scholes model, while for $\beta=-1$ this is the Bachelier, or time-
dependent Ornstein-Uhlenbeck (OU) model.. We assume that all parameters of the
model are known either as some continuous functions of time $t\in[0,\infty)$,
or as a discrete set of $N$ values for some moments $t_{i},i=1,\ldots,N$.

It is known that for the process $F_{t}$ to have a unique solution, an
explicit boundary condition has to be set for $-1<\beta<-1/2$ which is an
absorbing boundary condition for the price process to be a martingale and
arbitrage-free. The same boundary condition also holds naturally for
$-1/2\leq\beta<0$. In this case the SABR model has a probability mass at the
origin for $-1<\beta<0$ similar to that of the CEV model. In the $\sigma$
space it is relatively easy to show that for our SV model the boundary
$\sigma_{t}=0$ is an attainable regular boundary by Feller’s classification,
(Lipton, 2001). Therefore, similar to (Gorovoi and Linetsky, 2004) we always
make regular boundaries instantaneously reflecting, and include regular
reflecting boundaries into the state space. We also assume that infinite
boundaries are unattainable.

To attack the problem of pricing barrier options under this model we develop a
new method which is an extension of the Generalized integral transform (GIT)
method proposed in a series of the authors’ papers (some in cooperation with
Peter Carr and Alex Lipton) and covered in detail in a recent book (Itkin et
al., 2021). This method was originally developed to solve similar problems for
various time dependent one-factor models with time-dependent barriers and
provide semi-analytical prices of single and double barrier options. Then in
(Carr et al., 2021) for the first time this approach was extended to
stochastic volatility models. The authors developed the GIT method for pricing
barrier options in the time-dependent Heston model where the option price is
expressed in a semi-analytical form as a two-dimensional integral. This
integral depends on yet unknown function $\Psi(t,v),\ v=\sigma^{2}$ which is
the gradient of the solution at the moving boundary $F=L(t)$ and solves a
linear mixed Volterra-Fredholm (LMVF) equation of the second kind also derived
in that paper. Numerical examples illustrate high speed and accuracy of the
method as compared with the finite-difference approach. In this paper we
extend this idea and apply it to the time-dependent SABR model with mean-
reversion.

However, here we do it in a slightly different way. We represent the solution
of our problem as a weighted sum of the Bessel functions $J_{|\nu|}(x)$, so
actually this representation is a Fourier-Bessel series. The corresponding
weights solve a LMVF equation of the second kind which is also derived in the
paper. Thus, instead of the LMVF for $\Psi(t,v)$ we derive a LMVF equation for
the transform image ${\bar{u}}$. Once it is solved, the weights of the
Fourier-Bessel expansion become explicitly known, and so the solution is
expressed in closed form with no further integration. Also, in this way we are
able to keep time dependence of the model coefficients with no approximation,
while if we use the previous version of the GIT method, for this problem it
should be much harder if ever possible. Thus, this version of the GIT method
is another new result of this paper.

As shown in Section 5 (see also (Itkin et al., 2021)), our methods can be used
for any sort of barrier options. In this paper as an example we consider just
an Up-and-Out barrier Call option with $F_{t}\in[0,H(t)],\ H(t)>0$ is the
upper barrier. Once $F_{t}$ hits the barrier, the contract is terminated and
the option expires worthless with no rebate paid either at hit or at maturity,
i.e.

$C(t,H(t),\sigma)=0,$ (2)

where $C(t,F,\sigma)$ is the option price. This assumption is not restrictive
and can be relaxed, (Itkin and Muravey, 2021).

Note, that the barrier $H(t)$ is defined as a level of the forward price.
However, since it is a function of time, we can set the barrier for the spot
price $S_{t}$, so $S_{t}\in[H_{S}(t),\infty$ as well. In this case we have

$H(t)=H_{S}(t)\exp\left[\int_{0}^{t}(r(k)-q(k))dk\right],$

where $r(t),q(t)$ are the deterministic interest rate and continuous dividend.
In a similar way the model can be transformed to a similar model with respect
to the spot price $S_{t}$ by introducing a maturity dependent strike
$K_{S}=K\exp\left[\int_{0}^{T}(r(k)-q(k))dk\right]$.

At the other end of the domain we have $S=0$ and the Call price vanishes

$C(t,0,\sigma)=0.$ (3)

If the process $F_{t}$ survives till $t=T$, the option holder gets the Call
option payoff

$C(T,F,\sigma)=(F-K)^{+}.$ (4)

The Eq. (4) is the terminal condition for our problem. We also assume that
$H(T)\geq K$.

The rest of the paper is organized as follows. In Section 1 we consider the
PDE for the price of an Up-and-Out barrier Call option and solve it assuming
$\rho(t)=0$ by using the GIT method similar to that developed in (Carr et al.,
2020, 2021). In Section 2 we discuss how the LMVF equation derived in Section
1 can be solved numerically. In doing so we use the Radial Basis Functions
(RBF) method and show that using Gaussian RBFs makes the problem tractable by
reducing 2D integrals in the LVMF equation to the 1D ones. Section 3 extends
our method to the case when $\rho(t)\neq 0$. In Section 4 we describe some
numerical experiments where the barrier option prices obtained by using our
method and a FD approach are compared. We show that the speed and accuracy of
our method are comparable with those of the finite-difference approach at
small maturities and outperform them at high maturities even by using a
simplistic implementation of the RBF method.

We show that our method outperforms the FD one in both accuracy and speed.
Section 5 concludes.

##  1 The pricing PDE and its solution

The Feynman-Kac theorem, (Shreve, 1992), implies that under the risk neutral
measure the Call option price $C(t,F,\sigma)$ solves the partial differential
equation (PDE)

$\displaystyle\frac{\partial C}{\partial t}$
$\displaystyle+\frac{\displaystyle 1}{\displaystyle
2}\sigma^{2}F^{2(\beta+1)}\frac{\partial^{2}C}{\partial
F^{2}}+\frac{\displaystyle 1}{\displaystyle
2}\gamma^{2}(t)\sigma^{2}\frac{\partial^{2}C}{\partial\sigma^{2}}-\kappa(t)\sigma\frac{\partial
C}{\partial\sigma}+\rho(t)\gamma(t)\sigma^{2}F^{\beta+1}\frac{\partial^{2}C}{\partial
F\partial\sigma}=r(t)C,$ (5)

subject to the terminal condition in Eq. (4) and the boundary conditions in
Eq. (3). Eq. (2). Similar to (Carr et al., 2020), by making a change of
variables

$F=\left(-x\beta\right)^{-1/\beta},\qquad C(t,F,\sigma)\to
u(t,x,\sigma)e^{\int_{T}^{t}r(k)dk},$ (6)

we reduce this PDE to the form

$\displaystyle\frac{\partial u}{\partial t}$
$\displaystyle+\sigma^{2}\left[\frac{\displaystyle 1}{\displaystyle
2}\frac{\partial^{2}u}{\partial x^{2}}+\frac{b}{x}\frac{\partial u}{\partial
x}+\frac{\displaystyle 1}{\displaystyle
2}\gamma^{2}(t)\frac{\partial^{2}u}{\partial\sigma^{2}}+\rho(t)\gamma(t)\frac{\partial^{2}u}{\partial
x\partial\sigma}\right]-\kappa(t)\sigma\frac{\partial
u}{\partial\sigma}=0,\quad b=\frac{\beta+1}{2\beta}.$ (7)

When $\sigma$ is not stochastic Eq. (7) is the PDE associated with the one-
dimensional Bessel process, (Revuz and Yor, 1999)

$dX_{t}=dW_{t}+\frac{b}{X_{t}}dt.$ (8)

As mentioned in (Carr et al., 2020), to set the boundary and terminal
conditions in the new variables, we must distinguish two cases, which are
determined by the sign of $\beta$. If $-1<\beta<0$, the variable $x$ is
defined at $x\in[0,y(t)]$, with

$y(t)=-\frac{1}{\beta}H^{-\beta}(t)>0.$ (9)

Therefore, the boundary conditions now read

$u(t,0,\sigma)=u(t,y(t),\sigma)=0,$ (10)

and the terminal condition in Eq. (4) takes the form

$u(T,x,\sigma)=\left[\left(-\beta x\right)^{-1/\beta}-K\right]^{+},$ (11)

However, if $0<\beta<1$, the left boundary goes to $-\infty$. Therefore, in
this case it is convenient to redefine $x\to\bar{x}=-x$. Then $\bar{x}$ is
defined at $\bar{x}\in[\bar{y}(t),\infty)$ where

$\bar{y}(t)=\frac{1}{\beta}H^{-\beta}(t)>0.$ (12)

Accordingly, the terminal condition transforms to

$u(T,\bar{x},\sigma)=\left[\left(\beta\bar{x}\right)^{-1/\beta}-K\right]^{+},$
(13)

and the boundary conditions read

$u(t,\bar{x},\sigma)\Big{|}_{\bar{x}\to\infty}=u(t,\bar{y}(t),\sigma)=0.$ (14)

In other words, an Up-and-Out Call options written on the underlying process
$x$ in new variables behaves like a Down-and-Out Put option written on the
scaled underlying process $\bar{x}$ as in Eq. (12). It can also be checked
that when $0<\beta<1$ the PDE in Eq. (7) remains same but now in the $\bar{x}$
variable.

###  1.1 The GIT method for $\rho=0$

In this section we consider only the uncorrelated case when $\rho=0$. For the
sake of concreteness let us assume that $-1<\beta<0$. The other case
$0<\beta<1$ can be treated in a similar way, (Carr et al., 2020). Following
the approach of (Carr et al., 2020) for the CEV process with $x\in[0,y(t)]$ we
introduce the following integral transform

${\bar{u}}(t,\sigma,p)=\int_{0}^{y(t)}u(t,x,\sigma)x^{\nu+1}{J_{|\nu|}}(xp)dx,$
(15)

where $p=a+i\omega$ is a complex number with $-\pi/4<\arg(\sqrt{p})<\pi/4$ and
$\nu=1/(2\beta)<0$, since $\beta<0$. Multiplying both parts of Eq. (7) by
$x^{\nu+1}{J_{|\nu|}}(xp)$ and integrating on $x$ from $0$ to $y(t)$ we obtain

$\displaystyle\frac{\partial\bar{u}}{\partial t}$
$\displaystyle-y^{\prime}(t)[y(t)]^{\nu+1}{J_{|\nu|}}(y(t)p)u(t,y(t),\sigma)+\sigma^{2}\left[J_{1}+J_{2}+\frac{1}{2}\gamma^{2}(t)\frac{\partial^{2}\bar{u}}{\partial\sigma^{2}}\right]-\kappa(t)\sigma\frac{\partial\bar{u}}{\partial\sigma}=0.$
(16)

Here

$\displaystyle 2J_{1}$
$\displaystyle=\int_{0}^{y(t)}x^{\nu+1}{J_{|\nu|}}(xp)\frac{\partial^{2}u}{\partial
x^{2}}dx=x^{\nu+1}{J_{|\nu|}}(xp)\frac{\partial u}{\partial
x}\Bigg{|}_{0}^{y(t)}-u(t,x,\sigma)\frac{\partial}{\partial
x}\left(x^{\nu+1}{J_{|\nu|}}(xp)\right)\Bigg{|}_{0}^{y(t)}$ (17)
$\displaystyle+\int_{0}^{y(t)}u(t,x,\sigma)\frac{\partial^{2}}{\partial
x^{2}}\left(x^{\nu+1}{J_{|\nu|}}(xp)\right)dx,$ $\displaystyle J_{2}$
$\displaystyle=\int_{0}^{y(t)}x^{\nu+1}{J_{|\nu|}}(xp)\frac{b}{x}\frac{\partial
u}{\partial
x}dx=b\,x^{\nu}{J_{|\nu|}}(xp)u(t,x,\sigma)\Bigg{|}_{0}^{y(t)}-b\int_{0}^{y(t)}u(t,x,\sigma)\frac{\partial}{\partial
x}\left(x^{\nu}{J_{|\nu|}}(xp)\right)dx,$

With allowance for the boundary conditions Eq. (10), the sum $J_{1}+J_{2}$ can
be expressed as

$\displaystyle J_{1}+J_{2}=y^{\nu+1}(t){J_{|\nu|}}(y(t)p)\Psi(t,\sigma)$
$\displaystyle+\frac{1}{2}\int_{0}^{y(t)}u(t,x,\sigma)x\left[\frac{1-2\nu}{x}\frac{\partial}{\partial
x}\left(x^{\nu}{J_{|\nu|}}(xp)\right)+\frac{\partial^{2}}{\partial
x^{2}}\left(x^{\nu}{J_{|\nu|}}(xp)\right)\right]dx,$ $\displaystyle\Psi(t)$
$\displaystyle=\frac{1}{2}\frac{\partial u}{\partial x}\Big{|}_{x=y(t)},$ (18)

where we took into account that $b=\nu+1/2$.

Since the function $X(x)=x^{\nu}{J_{|\nu|}}(xp)$ solves the ordinary
differential equation (ODE), (Bateman and Erdélyi, 1953)

$\frac{d^{2}X}{dx^{2}}+\frac{1-2\nu}{x}\frac{dX}{dx}+p^{2}X=0,$ (19)

we can re-write Eq. (16) in the form

$\displaystyle\frac{\partial\bar{u}}{\partial t}$
$\displaystyle+\sigma^{2}\left[\frac{1}{2}\gamma^{2}(t)\frac{\partial^{2}\bar{u}}{\partial\sigma^{2}}-\frac{1}{2}p^{2}\bar{u}\right]-\kappa(t)\sigma\frac{\partial\bar{u}}{\partial\sigma}+\sigma^{2}y^{\nu+1}(t){J_{|\nu|}}(y(t)p)\Psi(t,\sigma)=0.$
(20)

The terminal condition for Eq. (20) follows from the definitions in Eq. (15)
and Eq. (11)

$\displaystyle{\bar{u}}(T,p)=-\frac{1}{p}\Bigg{\\{}$
$\displaystyle(-\beta)^{-1/\beta}\left[a^{1-\nu}J_{1-\nu}(ap)-y^{1-\nu}(T)J_{1-\nu}(py(T))\right]$
(21)
$\displaystyle-K\left[y^{1+\nu}(T)J_{-1-\nu}(py(T))-a^{1+\nu}J_{-1-\nu}(ap)\right]\Bigg{\\}},\qquad
a=-\frac{K^{-\beta}}{\beta}.$

Note, that since the option payoff doesn’t depend on $\sigma$, the terminal
condition in Eq. (21) either doesn’t depend on it.

The Eq. (20) is an inhomogeneous PDE in variables $(t,\sigma)$. In contrast to
the analogous PDE for the Heston model described and solved in (Carr et al.,
2021), this equation is not affine.

###  1.2 Solution of Eq. (20)

Our approach to solving Eq. (20) consists in few steps.

#### Step 1.

We begin with making a change of variables

$\displaystyle{\bar{u}}(t,\sigma,p)$ $\displaystyle=w(\tau,z,p),\quad
z(t,\sigma)=\log\sigma+g(t),\quad\tau(t)=\frac{1}{2}\int_{t}^{T}\gamma^{2}(k)dk,\quad
g(t)=\int_{T}^{t}\kappa(k)dk\,-\tau(t).$ (22)

In new variables Eq. (20) takes the form

$\displaystyle\frac{\partial w}{\partial\tau}$
$\displaystyle=\frac{\partial^{2}w}{\partial
z^{2}}+\eta(t,p)e^{2z}w+\frac{2}{\gamma^{2}(t)}e^{2(z-g(t))}y^{\nu+1}(t){J_{|\nu|}}(y(t)p)\Psi(t,z),$
(23) $\displaystyle t$
$\displaystyle=t(\tau),\quad\eta(t,p)=-\frac{p^{2}}{\gamma^{2}(t)}e^{-2g(t)}.$

Surprisingly, the homogeneous version of this PDE has been already considered
in (Itkin et al., 2020, 2021)444In (Itkin et al., 2020) we used the quadratic
drift, so instead of the term $n_{1}(t)z$ it gave rise to $n_{1}(t)e^{z}$.. We
have mentioned there that by the change of variables $z\to-z$ and
$\tau\to-\mathrm{i}\mkern 1.0mu\tau$, this PDE transforms into the time-
dependent Schrödinger equation with the unsteady Morse potential.

By using the Duhamel’s principle, one can deduce that the function $w$ solves
the LMVF equation of the second kind, (Itkin et al., 2020)

$\displaystyle w(\tau,z,p)$
$\displaystyle=\frac{1}{2\sqrt{\pi}}\Bigg{\\{}\int_{-\infty}^{\infty}\frac{w(0,\xi,p)}{\sqrt{\tau}}e^{-\frac{(z-\xi)^{2}}{4\tau}}d\xi+\int_{-\infty}^{\infty}d\xi\int_{0}^{\tau}\frac{w(k,\xi,p)}{\sqrt{\tau-k}}e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}}\eta(t(k),p)e^{2\xi}dk$
$\displaystyle+2\int_{-\infty}^{\infty}d\xi\int_{0}^{\tau}\frac{e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}}}{\gamma^{2}(k)\sqrt{\tau-k}}e^{2(\xi-g(k))}y^{\nu+1}(\tau){J_{|\nu|}}(y(k)p)\Psi(t(k),\xi)dk\Bigg{\\}},$
(24)

where $w(0,z,p)$ is the initial condition. Since
${\bar{u}}(T,\sigma,p)={\bar{u}}(T,p)$ doesn’t depend on $\sigma$ we must have

$w(0,z,p)={\bar{u}}(T,p).$ (25)

#### Step 2.

The inverse transform for Eq. (15) has been constructed in (Carr et al.,
2020). The solution is represented in the form

$u(t,x,\sigma)=x^{-\nu}\sum_{n=1}^{\infty}\alpha_{n}(t,\sigma)J_{|\nu|}\left(\frac{\mu_{n}x}{y(t)}\right),$
(26)

where $\mu_{n}$ is an ordered sequence of the positive zeros of
${J_{|\nu|}}(\mu)$:

$J_{|\nu|}(\mu_{n})=J_{|\nu|}(\mu_{m})=0,\quad\mu_{n}>\mu_{m}>0,\quad n>m.$

By definition, this form automatically respects the vanishing boundary
conditions for $u(t,x,\sigma)$ at $x\in[0,y(t)]$. We assume that under some
mild conditions on coefficients $\alpha_{n}(t,\sigma)$ this series converges
absolutely and uniformly $\forall x\in[0,y(t)]$ for any $t>0$.

The functions ${J_{|\nu|}}(\mu_{n}x),\ n=1,\ldots$ form an orthogonal basis in
the space $C[0,1]$ with the scalar product being

$\langle{J_{|\nu|}}(\mu_{k}x),{J_{|\nu|}}(\mu_{l}x)\rangle=2\int_{0}^{1}\frac{x{J_{|\nu|}}(\mu_{k}x){J_{|\nu|}}(\mu_{l}x)dx}{J_{|\nu|+1}(\mu_{k})J_{|\nu|+1}(\mu_{l})}=\begin{cases}1,\quad
k=l,\\\ 0,\quad k\neq l.\end{cases}$ (27)

Therefore, applying the direct integral transform in Eq. (15) to both parts of
Eq. (26) we obtain

$\bar{u}(t,\sigma,p)=y^{2}(t)\sum_{n=1}^{\infty}\alpha_{n}(t,\sigma)\int_{0}^{1}{\hat{x}}J_{|\nu|}\left(\mu_{n}{\hat{x}}\right){J_{|\nu|}}(py(t){\hat{x}})d{\hat{x}},\qquad{\hat{x}}=x/y(t).$
(28)

With allowance for Eq. (27), this equation can be solved for each coefficient
$\alpha_{n}(t,\sigma)$ to yield

$\alpha_{n}(t,\sigma)=2\frac{\bar{u}\left(t,\sigma,\mu_{n}/y(t)\right)}{y^{2}(t)J_{|\nu|+1}^{2}(\mu_{n})}.$
(29)

#### Step 3.

Using the inverse transform in Eq. (26) and coefficients
$\alpha_{n}(t,\sigma)$ found in Eq. (29), the solution $u(t,x,\sigma$ can be
represented as

$u(t,x,\sigma)=2x^{-\nu}\sum_{n=1}^{\infty}\frac{\bar{u}\left(t,\sigma,\mu_{n}/y(t)\right)}{y^{2}(t)J_{|\nu|+1}^{2}(\mu_{n})}J_{|\nu|}\left(\frac{\mu_{n}x}{y(t)}\right).$
(30)

Differentiating both parts of this representation by $x$, setting $x=y(t)$ and
using the identity $J_{|\nu|-1}(\mu_{n})=-J_{|\nu|+1}(\mu_{n})$, (Abramowitz
and Stegun, 1964), yields

$\Psi(t,\sigma)=-\frac{1}{y^{\nu+3}(t)}\sum_{n=1}^{\infty}\frac{\mu_{n}}{J_{|\nu|+1}(\mu_{n})}\bar{u}\left(t,\sigma,\mu_{n}/y(t)\right).$
(31)

Finally, substituting this expression into Eq. (1.2) we obtain a LMVF equation
of the second kind for $w(\tau,z,p)$

$\displaystyle w(\tau,z,p)$
$\displaystyle=\bar{u}(T,p)+\frac{1}{2\sqrt{\pi}}\int_{-\infty}^{\infty}d\xi\int_{0}^{\tau}\frac{w(k,\xi,p)}{\sqrt{\tau-k}}e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}}\eta(k,p)e^{2\xi}dk$
(32)
$\displaystyle-\frac{1}{\sqrt{\pi}}\sum_{n=1}^{\infty}\frac{\mu_{n}}{J_{|\nu|+1}(\mu_{n})}\int_{-\infty}^{\infty}d\xi\int_{0}^{\tau}\frac{y^{\nu+1}(\tau)w\left(k,\xi,\mu_{n}/y(k))\right)}{y^{\nu+3}(k)\gamma^{2}(k)\sqrt{\tau-k}}{J_{|\nu|}}(y(k)p)e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}+2(\xi-g(k))}dk.$

This equation has to be solved for each $p=p_{i}\equiv\mu_{i}/y(\tau),\
n=1,2,\ldots$. Once this is done, the solution of the whole problem is given
be Eq. (22). In particular, we are interested in
$\tau=\tau(0)=\frac{1}{2}\int_{0}^{T}\gamma^{2}(k)dk$.

It is easy to see that if $y(t)=y-const$ the last term in Eq. (32) vanishes,
and thus for each $w(\tau,z,\mu_{i}/y(\tau))$ we obtain an independent LMVF
equation

$\displaystyle w(\tau,z,\mu_{i}/y)$
$\displaystyle=\bar{u}(T,\mu_{i}/y)-\frac{\mu^{2}_{i}}{2\sqrt{\pi}y^{2}}\int_{-\infty}^{\infty}d\xi\int_{0}^{\tau}\frac{w(k,\xi,\mu_{i}/y)}{\gamma^{2}(k)\sqrt{\tau-k}}e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}+2(\xi-g(k))}dk,\quad
i=1,2,\ldots.$ (33)

Once they are solved, the solution of the whole problem is given be Eq. (22).

Looking at Eq. (32) one can discover that actually here we developed another
version of the GIT method which, to the best of the author’s knowledge, is not
known in the literature. Indeed, the first versions of the GIT method have
been proposed in (Kartashov, 2001) (see also reference therein), and as
applied to mathematical finance in (Carr and Itkin, 2021). Then, significant
development of this method for various problems and domains was done in a
series of the author’s papers in collaboration with Peter Carr and Alex
Lipton, see (Itkin et al., 2021). In all cases the core of the method is to
obtain a Volterra equation of the second kind for the gradient $\Psi(t)$ (in
the one-dimensional case), or the LMVF equation for $\Psi(t,\sigma)$ (in the
two dimensional case discussed in (Carr et al., 2021)). Once it is solved, the
solution is expressed in a semi-analytical form via an integral of $\Psi$.

Here we do it in a slightly different way. Instead of the LMVF equation for
$\Psi$ we derive a LMVF equation for the image ${\bar{u}}$. Once it is solved,
the solution is expressed in closed form (with no further integration). Also,
in this way we are able to keep time dependence of the coefficients
$\kappa(t),\gamma(t)$, while if we use the previous version of the GIT method
for solving this problem under consideration, this should be much harder if
ever possible.

###  1.3 Analysis of convergence

Recall, that $\nu=1/(2\beta)$ and we consider the case $-1<\beta<0$. A simple
analysis shows that the ratio $\frac{\mu_{n}}{J_{|\nu|+1}(\mu_{n})}$ rapidly
grows with the increase of $n$. To illustrate, in Fig. 1 the inverse ratio
$R_{n}(\nu)=\frac{J_{|\nu|+1}(\mu_{n})}{\mu_{n}}$ is depicted as a function of
$n$ for $\beta=-0.1,-0.4,-0.9$.

Figure 1: The ratio $R_{n}(\nu)=\frac{J_{|\nu|+1}(\mu_{n})}{\mu_{n}}$ as a
function of $n$ for $\beta=-0.1,-0.4,-0.9$.

To compensate this growth of the coefficients the solution
$w(\tau,z,\mu_{n}/y(\tau))$ should decrease faster than the coefficients to
converge. However, we don’t know in advance the speed of this decrease.
Therefore, numerical solution of Eq. (32) could be unstable. To get rid of
this effect we make a change of the dependent variable
$w(\tau,z,\mu_{n}/y(\tau))\mapsto{\bar{w}}(\tau,z,\mu_{n}/y(\tau))=w(\tau,z,\mu_{n}/y(\tau))/R_{n}(\nu)$.
Then Eq. (32) transforms to

$\displaystyle R_{i}(\nu){\bar{w}}(\tau,z,\mu_{i}/y(\tau))$
$\displaystyle=w(0,\mu_{i}/y(\tau))-\frac{\mu_{i}^{2}R_{i}(\nu)}{2\sqrt{\pi}}\int_{-\infty}^{\infty}d\xi\int_{0}^{\tau}\frac{{\bar{w}}(k,\xi,\mu_{i}/y(\tau))}{y^{2}(\tau)\gamma^{2}(k)\sqrt{\tau-k}}e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}+2(\xi-g(k))}dk$
(34)
$\displaystyle-\frac{1}{\sqrt{\pi}}\sum_{n=1}^{\infty}\int_{-\infty}^{\infty}d\xi\int_{0}^{\tau}\frac{y^{\nu+1}(\tau){\bar{w}}\left(k,\xi,\mu_{n}/y(k)\right)}{\gamma^{2}(k)y^{\nu+3}(k)\sqrt{\tau-k}}{J_{|\nu|}}\left(\mu_{i}\frac{y(k)}{y(\tau)}\right)e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}+2(\xi-g(k))}dk.$

Using Eq. (25) and Eq. (21), one can find that the first term in the RHS of
Eq. (34) has the form

$w(0,\mu_{i}/y(0))\propto\frac{J_{1-\nu}(\mu_{i})}{\mu_{i}}.$ (35)

With the increase of $i$ the roots $\mu_{i}\to\infty$. Therefore, in this
limit the function in Eq. (35) vanishes. Also, since $\nu=1/(2\beta)<0$
functions $w(0,\mu_{i}/y(0))$ in this limit rapidly tend to zero. Same is true
for coefficients $R_{i}(\nu)$ in the LHS of Eq. (34). Thus, assuming
${\bar{w}}(k,\xi,\mu_{i}/y(k))$ is finite for all $i=1,2,\ldots$, at high $i$
Eq. (34) asymptotically takes the form

$\displaystyle 0$
$\displaystyle=-\frac{1}{2\sqrt{\pi}}\int_{-\infty}^{\infty}d\xi\int_{0}^{\tau}dk\frac{e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}+2(\xi-g(k))}}{y^{2}(\tau)\gamma^{2}(k)\sqrt{\tau-k}}\left[{\bar{w}}(k,\xi,\mu_{i}/y(k))+2\frac{y^{\nu+1}(k)}{y^{\nu+1}(\tau)}\frac{{J_{|\nu|}}\left(\mu_{i}\frac{y(k)}{y(\tau)}\right)}{J_{|\nu|+1}(\mu_{i})\mu_{i}}\sum_{n=1}^{\infty}{\bar{w}}\left(k,\xi,\mu_{n}/y(k)\right)\right].$
(36)

With the increase of $i$ the term $P_{i}(\nu)=J_{|\nu|+1}(\mu_{i})\mu_{i}$
slowly grows, so $|J_{|\nu|+1}(\mu_{i})\mu_{i}|\to\infty,\ i\to\infty$.
Therefore, this forces the solution of Eq. (36) to vanish, i.e.
${\bar{w}}(k,\xi,\mu_{i}/y(k))\to 0,\,i\to\infty$.

Due to this behavior of ${\bar{w}}(k,\xi,\mu_{i}/y(k))$ and $P_{i}(\nu)$ at
high $i$, the solution in Eq. (30)

$\displaystyle u(t,x,\sigma)$
$\displaystyle=2x^{-\nu}\sum_{n=1}^{\infty}\bar{u}\left(t,\sigma,\mu_{n}/y(t)\right)\frac{1}{y^{2}(t)J_{|\nu|+1}^{2}(\mu_{n})}J_{|\nu|}\left(\frac{\mu_{n}x}{y(t)}\right)$
(37)
$\displaystyle=\frac{2x^{-\nu}}{y^{2}(t)}\sum_{n=1}^{\infty}\frac{{\bar{w}}\left(t,\sigma,\mu_{n}/y(t)\right)}{J_{|\nu|+1}(\mu_{n})\mu_{n}}J_{|\nu|}\left(\frac{\mu_{n}x}{y(t)}\right)$

converges. To illustrate the rate of convergence, in Eq. (37) we drop the
dependence of ${\bar{w}}\left(t,\sigma,\mu_{n}/y(t)\right)$ on $n$ and plot
$Z_{\nu}(M)=\sum_{n=1}^{M}J_{|\nu|}\left(\frac{\mu_{n}x}{y(t)}\right)/[J_{|\nu|+1}(\mu_{n})\mu_{n}]$
as a function of $M$ for $\beta=-0.1,-0.9$., where we use $\eta=x/y(t),\
0\leq\eta\leq 1$ and $\eta=0.3$. It can be seen in Fig. 2 that this function
becomes almost flat at $M>60$.

Figure 2: The behavior of $Z_{\nu}(M)$ as a function of $M$ for
$\beta=-0.1,-0.9$ where $\eta=x/y(t),\ 0\leq\eta\leq 1$ and $\eta=0.3$.

This behavior can be made even more transparent if we take into account, that
at $n\to\infty$

$\displaystyle\mu_{n}$ $\displaystyle\sim
n\pi+\frac{\pi}{2}\left(|\nu|-\frac{1}{2}\right)-\frac{4|\nu|^{2}-1}{8\pi(n+1/2(|\nu|-1/2))}+\mathcal{O}\left(\frac{1}{n^{3}}\right),$
(38) $\displaystyle J_{|\nu|}(\mu_{n})$
$\displaystyle\sim\left(\frac{2}{\pi\mu_{n}}\right)^{\frac{1}{2}}\left(\cos\omega\sum_{k=0}^{\infty}(-1)^{k}\frac{a_{2k}(|\nu|)}{\mu_{n}^{2k}}-\sin\omega\sum_{k=0}^{\infty}(-1)^{k}\frac{a_{2k+1}(|\nu|)}{\mu_{n}^{2k+1}}\right),$
$\displaystyle a_{k}(|\nu|)$
$\displaystyle=\frac{\left(4|\nu|^{2}-1^{2}\right)\left(4|\nu|^{2}-3^{2}\right)\cdots\left(4|\nu|^{2}-(2k-1)^{2}\right)}{k!8^{k}},\quad\omega=\mu_{n}-\frac{1}{2}\pi\left(|\nu|+\frac{1}{2}\right),$

see McMahon’s expansion in (Watson, 1966) and (Abramowitz and Stegun, 1964).
Taking the leading terms yields

$\displaystyle\frac{J_{|\nu|}\left(\frac{\mu_{n}x}{y(t)}\right)}{J_{|\nu|+1}(\mu_{n})\mu_{n}}$
$\displaystyle=\frac{(-1)^{1-n}\cos\left[\frac{1}{4}\pi(-2|\nu|+x(2|\nu|-4n+1)+1)\right]}{\pi
n\sqrt{x}}+O(n^{-2})$ (39) $\displaystyle=-\frac{1}{2\pi
n\sqrt{x}}\left\\{\cos\left[n\pi(1+x)-q\right]+\cos\left[n\pi(1-x)+q\right]\right\\}+O(n^{-2}),$
$\displaystyle q$ $\displaystyle=\frac{1}{4}\pi(2(x-1)|\nu|+x+1).$

Hence, at $n\to\infty$ the function $Z_{\nu}(M)$ behaves similar to the sum
$I_{1}+I_{2}$ of two cosine integrals

$I_{1}=\int_{n_{0}}^{\infty}\frac{\cos(n\pi(1+x)-q)}{n}\,dn,\qquad
I_{2}=\int_{n_{0}}^{\infty}\frac{\cos(n\pi(1-x)+q)}{n}\,dn,\quad x\neq 1,$
(40)

i.e. converges.

Since the function ${\bar{w}}\left(t,\sigma,\mu_{n}/y(t)\right)$ also
decreases with $n$ growing, the solution in Eq. (37) converges by the
Lebesgue’s dominated convergence theorem, and even faster than this is shown
in Fig. 2.

Recall, that the representation of the solution $u(t,x,\sigma)$ is given in
Eq. (30). One can observe that this semi-analytical solution can be treated as
the expansion of $u(t,x,\sigma)$ into Fourier-Bessel series with coefficients
proportional to function $\bar{u}\left(t,\sigma,\mu_{n}/y(t)\right)$. Each
such a function solves either the LMVF equation Eq. (33) (if $y(t)=y$ \-
const), or a system of LMVF equations in Eq. (34) (in a general case).

Alternatively, using the representation, (Watson, 1966)

$2\sum_{n=1}^{\infty}\frac{J_{|\nu|}\left(\frac{\mu_{n}x}{y(t)}\right)}{J_{|\nu|+1}(\mu_{n})\mu_{n}}=\left(\frac{x}{y(t)}\right)^{-\nu},\quad\nu<0,\
x\neq y(t),$ (41)

the solution in Eq. (37) can be bounded from above as

$\displaystyle u(t,x,\sigma)$
$\displaystyle=\frac{2x^{-\nu}}{y^{2}(t)}\sum_{n=1}^{\infty}\frac{{\bar{w}}\left(t,\sigma,\mu_{n}/y(t)\right)}{J_{|\nu|+1}(\mu_{n})\mu_{n}}J_{|\nu|}\left(\frac{\mu_{n}x}{y(t)}\right)\leq\frac{2x^{-\nu}}{y^{2}(t)}\max_{n}|{\bar{w}}\left(t,\sigma,\mu_{n}/y(t)\right)|\sum_{n=1}^{\infty}\frac{J_{|\nu|}\left(\frac{\mu_{n}x}{y(t)}\right)}{J_{|\nu|+1}(\mu_{n})\mu_{n}}$
(42)
$\displaystyle=\frac{x^{-2\nu}}{y^{2-\nu}(t)}\max_{n}|{\bar{w}}\left(t,\sigma,\mu_{n}/y(t)\right)|,\qquad
x\neq y(t).$

Since ${\bar{w}}\left(t,\sigma,\mu_{n}/y(t)\right)$ is the decreasing function
of $n$, and thus finite, the RHS of Eq. (42) is finite. Therefore, by the
Lebesgue’s dominated convergence theorem the solution of our problem is
finite.

It is, however, well known that Fourier–Bessel series are imposed to a Gibbs
Phenomenon. Also, even for continuous functions they could converge very slow,
(Gray and Pinsky, 1992). In our case this depends on the value of $\beta$. As
an example, let us put $\tau=0,\ y(t)=y=const$, so from Eq. (33) we obtain an
immediate solution

$w(0,z,\mu_{i}/y)=\bar{u}(T,\mu_{i}/y).$ (43)

Subsisting this solution into Eq. (30) we expect to get the option price equal
to the option payoff. We truncate infinite series in Eq. (30) up to $N$ terms
and compute the option price by summation. Convergence of the series in Eq.
(30) is depicted in Fig. 3 as a function of $N$ for several $\beta$. Other
parameters of the test are: $F=60,K=55$.

(a)  (b)

Figure 3: The relative error in % of the Call option value at maturity as a
function of the number of terms $N$ in the Fourier-Bessel series in Eq. (30)
for (a) $\beta=-0.1$; (b) - $\beta=-0.7$.

Fortunately, the Abelian theorem, (Hardy, 1991), states that for any summation
method $L$ if $u=\\{u_{n}\\}$ is a convergent sequence, with limit $U$ at
$n\to\infty$, then $L(u)=U$. In our case let us define $L$ to be the
arithmetic means of the first $N$ terms of $u$ (the Cesaro method). One can
prove that if $u$ does converge to $U$, then so does the sequence $Lu$. The
Call option price obtained by using the Cesaro method is also presented in
Fig. 3. It can be seen that this method significantly improves the
convergence, especially at low $|\beta|\ll 1$.

##  2 Solution of the LMVF equation Eq. (32)

It is known that solving the PDE in Eq. (5) can be done by using traditional
numerical methods, for instance the FD approach. As compared with this method,
solving the LMVF equation Eq. (34) is a bit less straightforward despite this
approach can have some advantages. First, for doing so we can use the RBF
method. We will show below that by using Gaussian RBFs the integrals on $\xi$
in Eq. (34) can be taken in closed form. Thus, instead of double integrals,
this equation will contain only single ones.

However, first we need to decide how to deal with the last term in Eq. (34)
which is new and doesn’t appear for the Heston model which is discussed in
(Carr et al., 2021). The problem with this term lies in the fact that it
contains summation over all functions ${\bar{w}}(k,\xi,\mu_{i}/y(k)),\
i=1,\ldots,$ and, thus, doesn’t allow factorization of Eq. (34) into a system
of independent LMVF for each $i\in[1,\infty)$. Therefore, we suggest solving
Eq. (34) iteratively.

The main idea behind this approach is as follows. As was mentioned in Section
1.3, functions ${\bar{w}}(k,\xi,\mu_{i}/y(k))$ tend to zero at high $i$.
Therefore, we can truncate the series and solve Eq. (34) only for $i\in[1,M]$,
where $M$ is finite and, perhaps, not so large. To begin, let us assume that
$W(k,\xi)=\sum_{i=1}^{\infty}{\bar{w}}(k,\xi,\mu_{i}/y(k))=0$. Then Eq. (34)
factorizes into $M$ independent LMVF equations (for each $i=[1,M]$) which are
given in Eq. (33).

Suppose this equation is solved numerically and the solution is denoted as
${\bar{w}}_{1}(\tau,z,\mu_{i}/y(\tau))$. Then at the next iteration the last
term in Eq. (34) can be replaced with the term

$\displaystyle\Lambda_{1}(\tau,z)=$
$\displaystyle-\frac{1}{\sqrt{\pi}}\sum_{n=1}^{\infty}\int_{-\infty}^{\infty}d\xi\int_{0}^{\tau}\frac{y^{\nu+1}(\tau){\bar{w}}_{1}\left(k,\xi,\mu_{n}/y(k)\right)}{\gamma^{2}(k)y^{\nu+3}(k)\sqrt{\tau-k}}{J_{|\nu|}}\left(\mu_{i}\frac{y(k)}{y(\tau)}\right)e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}+2(\xi-g(k)])}dk.$
(44)

Hence, instead of Eq. (34), at the second iteration we have to solve $M$
independent LMVF (for each $i\in[1,M]$ of the form

$\displaystyle{\bar{w}}_{2}\left(\tau,z,\frac{\mu_{i}}{y(\tau)}\right)$
$\displaystyle-\frac{\mu_{i}^{2}}{2\sqrt{\pi}}\int_{-\infty}^{\infty}d\xi\int_{0}^{\tau}\frac{{\bar{w}}_{2}(k,\xi,\frac{\mu_{i}}{y(\tau)})}{y^{2}(\tau)\gamma^{2}(k)\sqrt{\tau-k}}e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}+2(\xi-g(k))}dk={\bar{u}}\left(T,\frac{\mu_{i}}{y(0)}\right)+\frac{\Lambda_{1}(\tau,z)}{R_{i}(\nu)}.$
(45)

Once this is done, this process can be continued until convergence to the
desired tolerance in reached at some iteration $j$, so the final solution is
given by $w_{j}(\tau,z,\mu_{i}/y(\tau))$. The complexity of this process is
$jMC_{LMVF}$ where $C_{LMVF}$ is the complexity of solving a single LMVF of
the type Eq. (45). Since all LMVF are independent, this process is naturally
parallel and can be run on a parallel architecture very efficiently.

The LMVF of the type Eq. (45) can be solved by using, for instance, the RBF
method which we describe in the next section.

###  2.1 The RBF method

For the description of the RBF method as applied to the solution of integral
equations we refer the reader to (Assari et al., 2019; Zhang et al., 2014) and
references therein. Here we present only a short exposition of the method.

A function $\Phi:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is called to be radial
if there exists a univariate function $\phi:[0,\infty)\rightarrow\mathbb{R}$
such that

$\Phi(\mathbf{x})=\phi(r)$ (46)

where $r=\|\mathbf{x}\|$ and $\|\cdot\|$ is some norm in $\mathbb{R}^{d}$. In
this paper we consider just the Euclidean norm. Let
$\chi=\left\\{\mathbf{x}_{1},\ldots,\mathbf{x}_{N}\right\\}$ be a set of
scattered points selected in the domain $\Omega\subset\mathbb{R}^{d}$. A
function $u(\mathbf{x})$ at an arbitrary point $\mathbf{x}\in\Omega$ can be
approximated by using the global radial function $\phi(\|\mathbf{x}\|)$ via a
linear combination

$u(\mathbf{x})\approx\mathcal{G}_{N}u(\mathbf{x})=\sum_{i=1}^{N}c_{i}\phi\left(\left\|\mathbf{x}-\mathbf{x}_{i}\right\|\right),\quad\mathbf{x}\in\Omega,$
(47)

where the coefficients $\left\\{c_{1},\ldots,c_{N}\right\\}$ are determined by
the interpolation conditions

$\mathcal{G}_{N}u\left(\mathbf{x}_{i}\right)=u\left(\mathbf{x}_{i}\right)=u_{i},\quad
i=1,\ldots,N.$ (48)

In the literature various choices of the RBFs exist. Among others, here by a
tractability argument we prefer to use the Gaussian RBF

$\Phi(\mathbf{x})=e^{-\varepsilon\|\mathbf{x}\|^{2}},$ (49)

where $\varepsilon>0$ is the shape parameter. This function is strictly
positive-definite in $\mathbb{R}^{d}$ and, therefore, the expansion in Eq.
(47) is non-singular.

The advantage of using the Gaussian RBFs in our particular case lies in the
fact that we can eliminate one numerical integration in Eq. (34). Suppose we
are at the $m$-th iteration of the process described in the previous section.
Let us use Eq. (47) for interpolation of ${\bar{w}}(\tau,z,\mu_{i}/y(\tau))$
with the RBF given in Eq. (49)

${\bar{w}}(\tau,z,\mu_{i}/y(\tau))=\sum_{l}^{N_{l}}\sum_{j}^{N_{j}}c^{(i)}_{(j,l),m}e^{-\varepsilon(\tau-\tau_{j})^{2}-\varepsilon(z-z_{l})^{2}}.$
(50)

Substituting this expression into Eq. (34), after some algebra we obtain

$\displaystyle{\bar{u}}$
$\displaystyle(T,\mu_{i}/y(\tau))=\sum_{l}^{N_{l}}\sum_{j}^{N_{j}}c^{(i)}_{(j,l)}\left[e^{-\varepsilon(\tau-\tau_{j})^{2}-\varepsilon(z-z_{l})^{2}}+\mu_{i}^{2}\int_{0}^{\tau}\frac{e^{-2g(k)-\varepsilon(k-\tau_{j})^{2}}}{y^{2}(\tau)\gamma^{2}(k)}I_{jl}(z,\tau,k)\,dk\right],$
(51) $\displaystyle I_{jl}(z,\tau,k)$
$\displaystyle=\frac{1}{2\sqrt{\pi(\tau-k)}}\int_{-\infty}^{\infty}d\xi
e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}-\varepsilon(\xi-
z_{l})^{2}+2\xi}=\frac{1}{\sqrt{h}}\exp\left[-\frac{\epsilon(z-z_{l})^{2}-2z-4(\tau-k)(2z_{l}\epsilon+1)}{h}\right],$
$\displaystyle h$ $\displaystyle=h(\tau,k)=1+4\varepsilon(\tau-k).$

Here the term $\Lambda_{m}(\tau,z)$ can be represented as

$\displaystyle\Lambda_{m}(\tau,z)=$
$\displaystyle-\frac{1}{\sqrt{\pi}}\sum_{n=1}^{\infty}\sum_{l}^{N_{l}}\sum_{j}^{N_{j}}c^{(n)}_{(j,l),m}\int_{0}^{\tau}dk\frac{y^{\nu+1}(\tau)e^{-\varepsilon(k-\tau_{j})^{2}-2g(k)}}{\gamma^{2}(k)y^{\nu+3}(k)\sqrt{\tau-k}}{J_{|\nu|}}\left(\mu_{i}\frac{y(k)}{y(\tau)}\right)\int_{-\infty}^{\infty}d\xi
e^{-\frac{(z-\xi)^{2}}{4(\tau-k)}+2\xi-\varepsilon(\xi-z_{l})^{2}}$
$\displaystyle=-2\sum_{n=1}^{\infty}\sum_{l}^{N_{l}}\sum_{j}^{N_{j}}c^{(n)}_{(j,l),m}\int_{0}^{\tau}dk\frac{y^{\nu+1}(\tau)e^{-\varepsilon(k-\tau_{j})^{2}}}{y^{\nu+3}(k)\gamma^{2}(k)}{J_{|\nu|}}\left(\mu_{i}\frac{y(k)}{y(\tau)}\right)I_{jl}(z,\tau,k).$
(52)

The remaining integrals in $k$ in the RHS of Eq. (51) in general cannot be
taken in closed form and should be discretized by using some quadrature rule.
This is because the model coefficients $\kappa(t),\gamma(t)$ has to be found
by calibrating the model to market data. Alternatively, if practitioners are
able to guess the functional form of these coefficients, e.g.,
$\kappa(t)=\kappa_{1}-\kappa_{2}e^{-\kappa_{3}t},\ \kappa_{i}=const,\
i=1,2,3$, then they need to calibrate just the constants $\kappa_{i},\ i=1,3$.
In this case the integral on $k$ in Eq. (51) can be taken explicitly for some
families of functions, for instance for linear and quadratic functions of
$\tau$.

The RBF methods belong to the class of meshfree methods. That means that no
regular grid in $\tau_{j},z_{l}$ is required to run it (in contrast, e.g. to
the FD method). Therefore, taking a 2D set of collocation nodes
$\\{\tau_{j(l)},z_{l}\\},\ l=1,...,\bar{l},\ j(l)=1,...,\bar{j}(l)$ one can
substitute them into Eq. (51) and get a system of linear equations for the
coefficients $c_{j,l}$. For instance, in case of the regular grid with $N_{l}$
nodes in $z_{l}$ and $N_{j}$ nodes in $\tau_{j}$ we have $N_{j}N_{l}$ unknown
coefficients which solve the system of $N_{j}N_{l}$ linear equations. The
matrix of this system is dense, and therefore complexity of solving this
system by using the direct solver is $O(N^{3}_{j}N^{3}_{l})$. Obviously, this
result is not satisfactory from computational point of view. Iterative solvers
can improve this especially when a suitable preconditioner can be constructed.

It is, however, well-known, e.g., (Fasshauer and McCourt, 2012), the global
Gaussian RBF method leads to a notoriously ill-conditioned interpolation
matrix whenever $\varepsilon$ is small and the set of basis functions in Eq.
(49) becomes numerically linearly dependent on $\mathbb{R}^{d}$. This leads to
severe numerical instabilities and limits the practical use of Gaussians —
even though one can approximate a function with the Gaussian kernel with
spectral approximation rates. On the other hand, small $\varepsilon$ provide
better accuracy, and so have to be considered as an option for an accurate
pricing. It is also known that if $\varepsilon$ is kept fixed, convergence
stagnates even if $N_{j},N_{l}$ grow, and if $N_{j},N_{l}$ are fixed, the
error blows up with the decrease of $\varepsilon$.

One of the ideas to "fix" all these problems is using a "better basis" for RBF
interpolation to obtain well-conditioned (and therefore numerically stable)
interpolation, see (Fasshauer and McCourt, 2012) and references therein. Among
others, let us mention the RBF-QR method, (Fornberg et al., 2011; Larsson et
al., 2013; Larsson and Fornberg, 2005). Another approach would be instead of
the global RBF method using its localized version. The main improvement comes
from the fact that when two collocation nodes lie far away from each other,
the exponential function in Eq. (51) becomes small. Therefore, only the nodes
from the close neighborhood of the given node $(\tau_{j},z_{l})$ should be
taken into account in the summation in Eq. (51) (for instance, ${\cal N}$
closest nodes). This is similar to the FD method where numerical approximation
of derivatives is provided using just the closest nodes. In more detail, the
localized RBF method as applied to solving Volterra integral equations, is
described, e.g. in (Assari et al., 2019) (see also references therein).

When the localized version of RBF is used, the system matrix becomes sparse.
For a regular grid it becomes block-sparse with the size of the block equal to
$O({\cal N})$ (from ${\cal N}$ at edges to $2{\cal N}-1$ in the middle of the
matrix). Thus, the complexity of solving this system drops down to $O({\cal
N}^{2}N_{j},N_{l})$.

Having this in mind, in this paper we, however, use just a simple version of
the global Gaussian RBF method. This is done for two reasons. First, here we
want to illustrate that the proposed method gives reasonable option prices why
don’t require an immediate very high accuracy. Second, even with this global
method the speed of computations is better than that of the FD method. Also,
using more sophisticated RBF method for solving the LMVF equation thus
improving the speed and accuracy is subject of a separate paper which will be
presented elsewhere.

Going back to the global RBF method, we can re-write the interpolation
conditions in Eq. (47) in a matrix form

$\|B+A\mu_{i}^{2}\||C|=|{\bar{u}}|+|\Lambda_{m}|.$ (53)

Here $B$ is the $N_{j}N_{l}\times N_{j}N_{l}$ matrix
$\|e^{-\varepsilon(\tau-\tau_{j})^{2}-\varepsilon(z-z_{l})^{2}}\|$ \- the
Gaussian kernel for $\tau\in[\tau_{1},\ldots,\tau_{N_{j}}],\ \tau_{1}=0,\
\tau_{N_{j}}=\tau(0)$, $z\in[z_{1},\ldots,z_{N_{l}}]$, $A$ is a similar matrix
$\|I_{j},l\|$, ${\bar{u}}|$ is the $N_{j}N_{l}$ column vector with first
$N_{j}$ components equal to ${\bar{u}}(T,\mu_{i}/y(\tau_{1}))$, next $N_{j}$
components equal to ${\bar{u}}(T,\mu_{i}/y(\tau_{2}))$, and so on,
$|\Lambda_{m}|$ is the $N_{j}N_{l}$ column vector, and $C$ is the unknown
$N_{j}N_{l}$ column vector of coefficients. The Eq. (53) is a linear system of
algebraic equation whose solution gives us the coefficients $C$.

It is important to underline that for every $i$, the matrices $A,B$ are same,
as well as the matrix $\|B+A\mu_{i}^{2}\|$ doesn’t change in between the
iterations on $m$. Therefore, using a good matrix preconditioner helps to
reduce numerical complexity of iterations.

##  3 Non-zero correlation

When $\rho(t)\neq 0$ we solve the PDE in Eq. (7) by constructing a series in
$\rho(t)$ around the point $\rho(t)=0$. Thus, if the absolute value of the
correlation is small for all $t\geq 0$, several terms in the series expansion
could provide a suitable approximation to the exact solution. This kind of
expansion has been already considered in the literature for some SV models
with constant coefficients (see (Antonelli and Scarlatti, 2009) and references
therein). It was shown there that this power series converges with positive
radius under some regularity conditions.

Let us represent the solution of Eq. (7) as a series

$u(t,x,\sigma)=\sum_{i=0}^{\infty}\rho^{i}(t)u_{i}(t,x,\sigma),$ (54)

where $\rho(t)$ is a small parameter. Using this representation in the zero-
order approximation on $\rho(t)$ we get from Eq. (7)

$\displaystyle\frac{\partial u_{0}}{\partial t}$
$\displaystyle+\sigma^{2}\left[\frac{\displaystyle 1}{\displaystyle
2}\frac{\partial^{2}u_{0}}{\partial x^{2}}+\frac{b}{x}\frac{\partial
u_{0}}{\partial x}+\frac{\displaystyle 1}{\displaystyle
2}\gamma^{2}(t)\frac{\partial^{2}u_{0}}{\partial\sigma^{2}}\right]-\kappa(t)\sigma\frac{\partial
u_{0}}{\partial\sigma}=0.$ (55)

Let us also request that the first order approximation $u_{0}(t,x,\sigma)$
obeys the terminal condition in Eq. (11) and the boundary conditions in Eq.
(10). Then, it is easy to see that such a problem coincides with that for the
uncorrelated case. In Section 1.1 we obtained the solution of this problem
which is given by Eq. (37).

In the first order approximation in $\rho(t)$ from Eq. (7) we obtain the
following PDE

$\displaystyle\frac{\partial u_{1}}{\partial t}$
$\displaystyle+\sigma^{2}\left[\frac{\displaystyle 1}{\displaystyle
2}\frac{\partial^{2}u_{1}}{\partial x^{2}}+\frac{b}{x}\frac{\partial
u_{1}}{\partial x}+\frac{\displaystyle 1}{\displaystyle
2}\gamma^{2}(t)\frac{\partial^{2}u_{1}}{\partial\sigma^{2}}\right]-\kappa(t)\sigma\frac{\partial
u_{1}}{\partial\sigma}=-\gamma(t)\sigma^{2}\frac{\partial^{2}u_{0}}{\partial
x\partial\sigma}.$ (56)

As we have already satisfied the terminal and boundary conditions when
constructing a zero-order approximation, this PDE should be solved subject to
homogeneous terminal and boundary conditions.

The Eq. (56) is an inhomogeneous PDE. Suppose that the Green’s function
$G_{1}(t,x|T,x_{0})=G_{1}(x,x_{0},T-t)$ of the homogeneous counterpart of Eq.
(56) is known. Then, by using the Duhamel’s principle, (Polyanin, 2002), the
solution of Eq. (56) reads

$u_{1}(t,x,\sigma)=\int_{t}^{T}dk\int_{0}^{y(k)}d\xi\int_{0}^{\infty}dz\,G_{1}(T-k,x,\sigma,\xi,z)\gamma(k)z^{2}\frac{\partial^{2}u_{0}(k,\xi,z)}{\partial\xi\partial
z}.$ (57)

It is easy to see that next approximations can be constructed in a similar way
to provide

$u_{j}(t,x,\sigma)=\int_{t}^{T}dk\int_{0}^{y(k)}d\xi\int_{0}^{\infty}dz\,G_{j}(T-k,x,\sigma,\xi,z)\gamma(k)z^{2}\frac{\partial^{2}u_{j-1}(k,\xi,z)}{\partial\xi\partial
z}.$ (58)

Finally, observe that the homogeneous PDEs for all $u_{j}(t,x,\sigma),j\geq 0$
are same, i.e. $G_{j}=G_{0},\ \forall j\geq 0$. That means that we need to
find the Green’s function only once. By definition, it solves the PDE in Eq.
(55) with the boundary conditions in Eq. (10) and the terminal condition
$G(T,x,\sigma,\xi,z)=\delta(x-\xi)$ where $\delta(x)$ is the Dirac Delta
function. Thus, to find $G_{0}$ we can use the solution in Eq. (37) by just
replacing the terminal (initial in the variable $\tau$) condition.

##  4 Numerical experiments

Here we aim to test the accuracy and speed of the proposed approach by
comparing it with some FD method. For running numerical tests an explicit form
of the model parameters $\kappa(t),\gamma(t),\rho(t),\sigma,\beta$ should be
specified. Since we wish to investigate the speed and accuracy of the method,
there is no need to calibrate the model to real market quotes.

We compare the barrier option prices obtained by using our method with those
obtained by solving Eq. (5) using the FD method described in detail in (Itkin,
2015). This method belongs to a family of so-called ADI (alternative direction
implicit) schemes and provides second order of approximation in all
dimensions. Discretization of the PDE is done using a non-uniform grid which
is compressed close to the initial forward price $F$ and the initial
instantaneous volatility $\sigma$. A typical FD grid is presented in Fig. 4.
The method starts with doing few Rannacher to provide better stability of the
solution. The method has been validated for various problems. For instance, in
(Carr et al., 2021) it was used to price an European vanilla Put in the Heston
model with constant parameters
$\kappa_{0},\theta_{0},\sigma_{0},\rho_{0},v_{0}$ since for this model the Put
price can be found by FFT. It was observed that the relative error of the FD
scheme in use is about 12 bps.

Also, for the Call option at $t=T$ there is a kink in the option price at the
upper boundary $H$. Therefore, for a better stability it is useful to solve
this problem for a covered Call ${\cal C}(t,F)=C(t,F)-F+K$. For the covered
Call the boundary conditions obviously change to ${\cal C}(t,0)=K,\ {\cal
C}(t,H)=K-H$ and the terminal condition - to ${\cal C}(T,F)=(F-K)^{+}+K-F$.
Thus, when the Call is In-The-Money at $t=T$, the payoff of the covered Call
vanishes. Still at $F=H,\ t=T$ there is a kink in the payoff of the covered
Call, but now of a much less amplitude. Also, the PDE for ${\cal C}(t,F)$
coincides with that for $C(t,0)$ but now has an extra source term $-r(F-K)$.
This term, however, doesn’t change the ADI scheme, and can be added to the
explicit steps of the methods.

Figure 4: A typical nonuniform grid with 76 nodes in $F$ and 79 nodes in
$\sigma$ for the FD scheme with $F_{0}=60,\ \sigma_{0}=0.5$ and $K=50$.

#### Example 1.

In this test we assume that the instantaneous volatility $\sigma$ is constant,
i.e. not stochastic. Also let $y(t)=y-const,\ r(t)=r-const$. In this case Eq.
(20) can be solved analytically which yields

${\bar{u}}(t,\sigma,p)={\bar{u}}(T,p)e^{-\frac{1}{2}p^{2}\sigma^{2}(T-t)}.$
(59)

Thus, the Up-and-Out barrier Call option price can be expressed analytically
as

$C(x,\sigma)=2e^{-rT}x^{-\nu}\sum_{n=1}^{\infty}{\bar{u}}(T,\mu_{n}/y)e^{-\frac{\mu^{2}_{n}}{2y^{2}}\sigma^{2}T}\frac{J_{|\nu|}\left(\frac{\mu_{n}x}{y}\right)}{J^{2}_{|\nu|+1}(\mu_{n})}.$
(60)

Note, that this expression can be also represented in a different form
(despite, perhaps, less practical). Indeed, taking into account the definition
of ${\bar{u}}(T,p)$ as the GIT of the option payoff and substituting it into
Eq. (60) we get

$\displaystyle C(x,\sigma)$
$\displaystyle=2e^{-rT}x^{1-2\nu}y\int_{0}^{1}\left[A\xi^{-2\nu}-K\right]^{+}\Theta_{|\nu|}(\sigma\sqrt{T}/y,\xi,\xi)d\xi,\qquad
A=\left(-\beta y\right)^{-1/\beta},$ (61)

where the function $\Theta_{|\nu|}(\varsigma,x_{1},x_{2})$ was introduced in
(Carr et al., 2020) as

$\Theta_{|\nu|}(\varsigma,x_{1},x_{2})=\sum_{n=1}^{\infty}e^{-\frac{\mu^{2}_{n}\varsigma}{2}}\frac{x_{1}^{\nu}J_{|\nu|}(\mu_{n}x_{1})}{J_{|\nu|+1}(\mu_{n})}\frac{x_{2}^{\nu}J_{|\nu|}(\mu_{n}x_{2})}{J_{|\nu|+1}(\mu_{n})}.$
(62)

This function (let us call it the Bessel Theta function) is a Bessel analog of
the Jacobi Theta function of the $\theta_{3}(\phi,x)$, (Mumford et al., 1983),
in a sense, that $\theta_{3}(\phi,x)$ is a periodic solution of the heat
equation, while $\Theta_{|\nu|}(\varsigma,x_{1},x_{2})$ is a periodic solution
of the Bessel equation. It can be checked that, if $\nu=1/2$, we have

$\displaystyle\Theta_{1/2}\left(\frac{\sqrt{\tau}}{y(\tau)},\frac{s}{y(\tau)},\frac{z}{y(\tau)}\right)=\frac{y(\tau)}{2}\left[\theta_{3}\left(e^{-\frac{\pi^{2}\theta^{2}}{8}},\frac{n\pi(s-z)}{4y(\tau)}\right)-\theta_{3}\left(e^{-\frac{\pi^{2}\theta^{2}}{8}},\frac{n\pi(s+z)}{4y(\tau)}\right)\right].$
(63)

Also, it can be checked that the Bessel Theta function at $\varsigma\to 0$
becomes the Dirac delta function. Therefore,
$\Theta_{|\nu|}(\varsigma,x_{1},x_{2})$ in Eq. (61) is a scaled density of the
underlying process at the domain $F\in[0,y]$, and Eq. (61) can be naturally
treated as a discounted expectation of the option payoff under risk-neutral
measure.

However, as we mentioned, representation Eq. (61) is less practical. That is
because the Jacobi theta functions are implemented in many numerical libraries
and programming languages, like MpMath in python, or Wolfram Mathematica or
Matlab, etc., while the Bessel Theta function is not. Therefore, manual
implementation is required at the moment.

To proceed with the analytical representation of the Call option price in Eq.
(60), we take parameters of the model as in Table 1 and look at convergence of
the option price as a function of the number of terms $N$ in Eq. (60). The
results are presented in Fig. 5. It can be seen that all series converge at
$N\approx 150$ for short maturities, and even faster for longer maturities.

$F$ | $\sigma$ | $T$ | $r$ | $H$ | $K$ | $\beta$
---|---|---|---|---|---|---
60 | 0.5 | 1/12, 0.25 | 0.02 | 80 | 55 | -0.1, -0.7
Table 1: Parameters of the model in Example 1.

Figure 5: The Up-and-Out barrier Call option price computed by using Eq. (60)
as a function of $N$ \- the number of terms in the series.

In Table 2 we present the comparison of results computed by using Eq. (60)
with $N=250$ with those obtained by using the FD method. As here the
instantaneous volatility is constant and so the model becomes one-dimensional,
we also provide the results obtained by using a 1D FD scheme described in
(Carr et al., 2020) and also elapsed times of all methods.

Method | T, yrs
---|---
| 0.042 | 0.083 | 0.250 | 0.500 | 1.000 | 2.000 | El. time, sec
Analytic | 5.1816 | 5.4908 | 5.0768 | 3.5365 | 1.8997 | 0.8333 | 0.007
FD 1D | 5.1811 | 5.4875 | 5.0731 | 3.5352 | 1.8994 | 0.8332 | 0.17
FD 2D | 5.1805 | 5.4813 | 5.0617 | 3.5299 | 1.8975 | 0.8326 | 1.60

Table 2: The Up-and-Out barrier Call option prices for $\sigma$ = const,
$\beta$ = -0.1, computed by various methods and the elapsed time in secs (for
one cell in the table).

It can be seen that all prices obtained by various methods are close to each
other while the analytic computation is obviously much faster.

Same results for $K=60$ are presented in Table 3.

Method | T, yrs
---|---
| 0.042 | 0.083 | 0.250 | 0.500 | 1.000 | 2.000 | El. time, sec
Analytic | 1.6203 | 2.2467 | 2.5756 | 1.8174 | 0.9565 | 0.4104 | 0.007
FD 1D | 1.6172 | 2.2451 | 2.5729 | 1.8164 | 0.9562 | 0.4103 | 0.17
FD 2D | 1.6182 | 2.2427 | 2.5647 | 1.8085 | 0.9476 | 0.4043 | 1.60

Table 3: The Up-and-Out barrier Call option prices for $\sigma$ = const,
$\beta$ = -0.1 and $K=60$ (ATM) computed by various methods and the elapsed
time in secs (for one cell in the table).

#### Example 2.

Here we consider the whole $\lambda$-SABR model, so the instantaneous
volatility is stochastic. Since we don’t calibrate the model, without loss of
generality we choose just some artificial time dependencies of the model
parameters, namely:

$\gamma(t)=\gamma_{1}e^{-\gamma_{2}t},\quad\rho(t)=0,\quad\kappa(t)=\kappa_{1}e^{-\kappa_{2}t},\quad
r(t)=r_{0},\quad H(t)=H.$ (64)

where $\gamma_{1},\gamma_{2},\kappa_{1},\kappa_{2}$ are constants. With these
definitions we can directly find

$\displaystyle\tau(t)$
$\displaystyle=\frac{\gamma_{1}^{2}}{4\gamma_{2}}\left(e^{-2\gamma_{2}t}-e^{-2\gamma_{2}T}\right),\quad
g(t)=\frac{\kappa_{1}}{\kappa_{2}}\left(e^{-\kappa_{2}T}-e^{-\kappa_{2}t}\right)-\tau(t),\quad\eta(t,p)=-\frac{p^{2}}{\gamma^{2}(t)}e^{-2g(t)}.$
(65)

Other parameters of the test are presented in Table 4.

$F$ | $\sigma$ | $\gamma_{1}$ | $\gamma_{2}$ | $\kappa_{1}$ | $\kappa_{2}$ | $H$ | $r_{0}$
---|---|---|---|---|---|---|---
60 | 0.5 | 0.5 | 0.3 | 1. | 0.2 | 80 | 0.02
Table 4: Parameters of the model in Example 2.

We run the test for a set of maturities $T\in[1/24,1/12,0.25,0.5,1,2]$ years
and strikes $K\in[45,50,55,60,65,70]$. For the RBF method as the collocation
points we use a uniform grid in $t\in[0,T]$ and
$z\in[z_{0}-z_{m},z_{0}+z_{m}]$ where $z_{0}=\log\sigma+g(0)$, $z_{m}$ = 0.5.
Close to ATM it is useful to add an extra (or even few) collocation point
behind $\tau(0)$ for better approximation. We take $N_{j}=10,N_{l}=20$,
$N=350$. The best values of $\varepsilon$ in all experiments are given in
Table 5. It can be seen that these values are almost independent of the
strikes and maturities, as well as of $\beta$. Only for the ATM options we
need to decrease $\varepsilon$555Since the payoff function has a kink ATM, it
is always a problem for any numerical method to address it in the numerical
solution. Various approached was developed in the literature, see a survey,
e.g., in (Itkin, 2017)..

$\mathbf{\beta}$ | K
---|---
| 45 | 50 | 55 | 60 | 65 | 70
-0.1 | 0.10 | 0.10 | 0.10 | 0.02 | 0.15 | 0.15
-0.7 | 0.15 | 0.15 | 0.15 | 0.02 | 0.15 | 0.10
Table 5: The values of $\varepsilon$ used in the numerical experiments.

When solving Eq. (51) we approximate the integral in time by using a Simpson
method with $N$ (in general, non-even) nodes distributed in $k\in[0,\tau]$.
For the FD method the time step is fixed and equal to $\min(0.01,T/50)$ year
to preserve the method’s accuracy in time at high $T$. Accordingly, since our
FD scheme is of the second order of approximation in time, and the Simpson
rule provides the fourth order of approximation, to preserve same accuracy
using the Simpson method we need to take $N_{j}\approx 10T$. This dictates our
choice of $N_{j}=10$ for all maturities (numerical experiments show that even
for $T=2$ the increase of $N_{j}$ doesn’t improve the results).

To solve Eq. (53) a standard `bicgstabl` iterative solver is used with no
preconditioner. The Up-and-Out barrier Call option prices computed in these
experiments using Matlab are presented in Tab. 6 and also in Fig. 6. Typical
elapsed times are also shown in Tab. 6. With $N=250$ the accuracy of the RBF
method slightly drops down while the elapsed time becomes about 1.1 sec.

$\beta=-0.1$ $\bf{T}$ | 0.038 | 0.083 | 0.25 | 0.5 | 1 | 2  | 0.038 | 0.083 | 0.25 | 0.5 | 1 | 2
---|---|---|---|---|---|---|---|---|---|---|---|---
$\bf{K}$ |  GIT | FD
70  | 0.0412 | 0.1151 | 0.2815 | 0.3452 | 0.2456 | 0.1997  | 0.0149 | 0.0992 | 0.2713 | 0.2635 | 0.2255 | 0.2005
65  | 0.3286 | 0.6547 | 0.9666 | 0.9584 | 0.6648 | 0.6847  | 0.2280 | 0.5799 | 1.0118 | 0.9352 | 0.7918 | 0.7027
60  | 1.7614 | 2.2653 | 2.6923 | 2.3719 | 1.8738 | 1.7136  | 1.5838 | 2.1574 | 2.5960 | 2.2809 | 1.9143 | 1.6972
55  | 5.1455 | 5.3522 | 5.2544 | 4.4732 | 3.6379 | 3.3626  | 5.1691 | 5.4385 | 5.2927 | 4.5130 | 3.7610 | 3.3325
50  | 9.7614 | 9.6848 | 9.1691 | 7.4481 | 6.1788 | 5.6118  | 9.9943 | 9.9678 | 9.0054 | 7.5949 | 6.3409 | 5.6325
45  | 14.6288 | 14.4850 | 13.2513 | 11.0707 | 9.2756 | 8.4581  | 14.9845 | 14.8887 | 13.3691 | 11.3469 | 9.5560 | 8.5304
Elapsed time  | 1.5 | 1.5 | 1.5 | 1.5 | 1.5 | 1.5  | 1.5 | 1.5 | 1.5 | 1.5 | 3.1 | 6.1
$\beta=-0.7$ $\bf{T}$ | 0.038 | 0.083 | 0.25 | 0.5 | 1 | 2  | 0.038 | 0.083 | 0.25 | 0.5 | 1 | 2
$\bf{K}$ |  GIT | FD
70  | 0.0038 | 0.0136 | 0.0131 | 0.0254 | 0.0126 | 0.0027  | 6.9416E-08 | 2.7755E-07 | 2.4938E-06 | 1.0001E-05 | 2.1420E-05 | 5.1659E-05
65  | 0.0130 | 0.0156 | 0.0237 | 0.0146 | 0.0271 | 0.0037  | 3.4708E-08 | 1.3877E-07 | 1.3181E-06 | 1.6641E-05 | 2.0137E-04 | 6.7151E-04
60  | 0.2432 | 0.2854 | 0.3486 | 0.4184 | 0.4538 | 0.4711  | 0.1258 | 0.1806 | 0.2968 | 0.3800 | 0.4494 | 0.4820
55  | 5.0148 | 4.9903 | 4.9794 | 4.9982 | 4.9096 | 4.7928  | 4.9958 | 4.9917 | 4.9751 | 4.9503 | 4.9011 | 4.8045
50  | 9.9764 | 9.9895 | 9.9332 | 9.8322 | 9.7938 | 9.6076  | 9.9917 | 9.9833 | 9.9501 | 9.9005 | 9.8020 | 9.6079
45  | 14.9749 | 14.9291 | 14.8793 | 14.8485 | 14.7052 | 14.4162  | 14.9875 | 14.9750 | 14.9252 | 14.8507 | 14.7030 | 14.4118
Elapsed time  | 1.5 | 1.5 | 1.5 | 1.5 | 1.5 | 1.5  | 1.6 | 1.4 | 1.4 | 1.4 | 2.7 | 5.7

Table 6: Comparison of Up-and-Out Call option prices for the $\lambda$-SABR
model obtained by various methods and elapsed time in secs.

Figure 6: Up-and-Out Call option prices for the $\lambda$-SABR model obtained
by the FD method in our experiment, with $\beta$ = -0.1.

Fig. 7 presents a typical Up-and-Out barrier Call option surface $C(F,\sigma)$
obtained for $\beta=-0.1,\ K=55,\ T=1/24$, and other model parameters as in
Table. 6.

Figure 7: Up-and-Out barrier Call option $C(F,\sigma)$ for $\beta=-0.1,\
K=55,\ T=1/24$, and other model parameters as in Table. 6.

In Table 7 the percentage error between the option prices obtained by our and
the FD method are presented as a function of $K,T$.

$\beta=-0.1$ | T
---|---
K | 0.0380 | 0.0830 | 0.2500 | 0.5000 | 1.0000 | 2.0000
70 | 63.74 | 13.78 | 3.61 | 23.66 | 8.20 | -0.41
65 | 30.63 | 11.42 | -4.67 | 2.42 | -19.09 | -2.63
60 | 10.08 | 4.76 | 3.58 | 3.84 | -2.16 | 0.96
55 | -0.46 | -1.61 | -0.73 | -0.89 | -3.38 | 0.89
50 | -2.39 | -2.92 | 1.78 | -1.97 | -2.62 | -0.37
45 | -2.43 | -2.79 | -0.89 | -2.50 | -3.02 | -0.85
$\beta=-0.7$ | T
K | 0.0380 | 0.0830 | 0.2500 | 0.5000 | 1.0000 | 2.0000
70 | 100.00 | 100.00 | 99.98 | 99.96 | 99.83 | 98.10
65 | 100.00 | 100.00 | 99.99 | 99.89 | 99.26 | 81.93
60 | 48.29 | 36.74 | 14.88 | 9.16 | 0.98 | -2.30
55 | 0.38 | -0.03 | 0.09 | 0.96 | 0.17 | -0.24
50 | -0.15 | 0.06 | -0.17 | -0.69 | -0.08 | 0.00
45 | -0.08 | -0.31 | -0.31 | -0.01 | 0.02 | 0.03
Table 7: The percentage error between the option prices obtained by our and
the FD method as a function of $K,T$.

Big relative error at $K=70,65$ doesn’t seem to confuse readers since the
option price itself is very small (several cents at $K=70$). However, the ATM
prices have bigger error. We attribute this error to both the simplistic
numerical approach we use to solve the LMVF equation, and to the accuracy of
the FD solver for the ATM prices. This can be seen from Tables 2,3 where the
difference between the ATM prices computed by our method with the exact
solution for ${\bar{u}}$ (when $\gamma(t)=\kappa(t)=0$) and those obtained by
using the FD methods also slightly increases as compared with the OTM prices.

The agreement between various methods is much better for the ITM prices. This
is because those are closer to the intrinsic value and so the values of the
integrals are small (and so are the errors). Also, the errors are bigger at
small maturities and lower at the high ones. A partial explanation of this
effect was provided in (Carr et al., 2020), where the authors compared the GIT
and HP (heat potential) methods as applied to pricing barrier options under
one-dimensional CIR and CEV models. It was observed that the GIT method
provides very accurate results at high maturities in contrast to the HP method
which is good at small $T$. This can be verified by looking at exponents under
the GIT solution integral (see, e.g., Eq. (59)) which are proportional to $T$.
Contrary, for short maturities the GIT method is slightly less accurate than
the HP method, as the exponents in the HP solution integral are inversely
proportional to $T$.

Another possible reason could be different properties of the RBF matrix at
various maturities. To check that we recomputed the ATM prices by using a
`minres` solver which is good when the matrix is not positive definite, but
symmetric, and we can still construct an orthogonal basis for the Krylov
subspace by three-term recurrence relations, (Paige and Saunders, 1975).
Indeed, we do see that due to small rounding errors some eigenvalues of our
RBF matrices are very small negative numbers, therefore, in principle, the RBF
matrix should be first regularized by using a standard procedure. However,
again, since this error depends on quality of the RBF interpolation it is
expected that modern methods which are more stable than the Gaussian global
RBF method would provide better accuracy. Indeed, for the Gaussian RBF method
a typical matrix in Eq. (53) has the condition number about $10^{19}$ and,
therefore, even iterative methods of solving this system of linear equations
produce bigger error than in the case of a well-conditioned matrix.

We also double the value of $N$ when using `minres` at $\beta=-0.7$ while
doing so for `bicgstabl` doesn’t change the computed prices. The results of
these tests for the ATM options are presented in Table 8. It can be seen that
using the `minres` solver significantly drops down the relative error as
compared with the `bicgstabl` solver at small and intermediate maturities,
while the later solver outperforms the former at high maturities. Also
`minres` is almost twice faster than `bicgstabl`, and so twice faster than the
FD solver at small and intermediate $T$.

$\beta$ / solver | T, C | T, % dif
---|---|---
| 0.0380 | 0.0830 | 0.2500 | 0.5000 | 1.0000 | 2.0000 | 0.0380 | 0.0830 | 0.2500 | 0.5000 | 1.0000 | 2.0000
-0.7/bicgstabl | 0.2432 | 0.2596 | 0.3634 | 0.4055 | 0.4094 | 0.3919 | 48.2942 | 30.4523 | 18.3371 | 6.2793 | -9.7604 | -22.9764
-0.7/minres | 0.1503 | 0.1869 | 0.2644 | 0.3429 | 0.3685 | 0.3331 | 16.3537 | 3.3986 | -12.2616 | -10.8419 | -21.9523 | -44.6957
-0.1/bicgstabl | 1.7614 | 2.2653 | 2.6923 | 2.3719 | 1.8738 | 1.7136 | 10.0814 | 4.7629 | 3.5769 | 3.8370 | -2.1570 | 0.9576
-0.1/minres | 1.6073 | 2.0037 | 2.2804 | 2.2098 | 2.1352 | 2.1290 | 1.4603 | -7.6687 | -13.8411 | -3.2192 | 10.3477 | 20.2828

Table 8: The Up-and-Out ATM Call option price computed by using two different
solvers for $\beta$ = -0.7, -0.1 (left panel) and the percentage difference
between these prices and those computed by using the FD method (right panel).

##  5 Discussion

In this paper we extend the GIT method proposed to get semi-analytical prices
of barrier options for various 1D time-dependent models with time-dependent
barriers, see (Itkin et al., 2021) and references therein, by applying it to
the so-called $\lambda$-SABR stochastic volatility model. In doing so we
develop a new modification of the GIT method which delivers solution of this
problem as a Fourier–Bessel series. Coefficients of this series solve a LMVF
equation of the second kind which is derived and discussed in the paper.

Let us compare this approach step-by-step with the numerical FD method.

  1. 1.

With the FD method we create a 3D grid in $x-\sigma-t$ space, so the
complexity of the solution is about $O(N_{t}N_{x}N_{\sigma})$, where $N_{i}$
is the corresponding number of nodes in the $i$-th direction. With the GIT
method we still some nodes in the $t$ and $\sigma$ space. Since integration in
the $t$ space can be done by using high order quadratures, the number of the
RBF nodes $n_{t}$ can be significantly reduced providing same accuracy. Thus,
in general $n_{t}\ll N_{t}$. In the $\sigma$ space the GIT method doesn’t
require numerical integration (it is done analytically), and thus, we
typically have $n_{\sigma}\ll N_{\sigma}$.

  2. 2.

The GIT method doesn’t need any grid in $x$. Instead, it needs to compute
multiple terms in the Fourier-Bessel series. Fortunately, all matrices in Eq.
(53) can be precomputed and reused. Therefore, solution of these systems for
various $\mu_{i},i=1,\ldots,$ is fast. Since the series converges, it can be
truncated to some fixed number of terms.

Also, all terms in this series are produced by solving the corresponding LMVF
equation which are independent of each other. Therefore, they can be solved in
parallel by using any parallel architecture.

  3. 3.

We compared pricing of Call barrier options by using the GIT and FD methods
strike-by-strike and maturity-by-maturity. In other words, we used a backward
approach. If one needs prices for multiple $K$ and $T$ a better way would be
to solve the forward equation for the density of the underlying process and
then integrate the solution with payoffs thus computing the discounted
expectation. It turns out that the proposed GIT method is a perfect fit for
doing so. Indeed, when solving the LMVF equation we can take $\tau$ to be
$\tau(\max(T))$. Since our method uses the RBF interpolation, it generates the
results for all $\tau\in[0,\tau(\max(T))]$ just in one sweep. Therefore,
multiple $T$ can be priced at once. For multiple strikes, assuming that
$\varepsilon$ is independent on $K$, the only remaining function of $K$ is
${\bar{u}}(T,p)$ which enters the RHS of Eq. (53). Therefore, it can be very
fast precomputed for multiple strikes (just a vector operation), and then we
can solve Eq. (53) for various RHS simultaneously. Indeed, it requires just
solving one system with multiple RHS sides that could be efficiently done with
the modern software. Thus, prices for all pairs $T,K$ can be obtained in one
sweep. Therefore, the relative performance of the forward FD and GIT methods
remains the same (the FD method also requires computing multiple expectations,
i.e. the integrals with the lower limit being equal to $K$).

  4. 4.

Based on the results obtained in the paper, performance of the GIT method is
almost twice better than that of the FD method at small maturities (by using
the `minres` solver), of the same order at intermediate maturities (by using
the `bicgstabl` solver), while the GIT method is faster at long maturities.
The accuracy of the GIT method combined together with the Gaussian RBF is good
but depends on the value of the shape parameter $\varepsilon$ (as it should
be) and a choice of the linear solver. We show that choosing an appropriate
solver gives rise to good agreement between our and FD benchmark results.

  5. 5.

To solve the LMVF equations we used the Gaussian RBF method. As mentioned in
Section 2.1, modern RBF approaches significantly improve both speed and
accuracy of the solution, and, therefore, it is expected that using these
methods will significantly improve both speed and accuracy. This is subject of
our future investigation which we will report elsewhere.

On a general note, solving a problem with time-dependent barriers by using a
FD method brings some technical challenges. In this case either the FD gird
has to be reconstructed at every time step to have the barrier as the last
node, or a fixed grid can be used but the boundary condition should be moved
along the grid, e.g., by interpolating the barrier onto the closest node. The
later decreases the accuracy of the method while the former deteriorates the
speed. In contrast, our semi-analytical method doesn’t face this problem at
all.

Another advantage of our method as compared with the FD method is that it is
semi-analytical. This means, that various properties of the solution can be
retrieved by analyzing its Fourier-Bessel series representation. For instance,
computation of Greeks can be done by differentiating this series by the
necessary parameter of the model. In particular, derivatives on $x$ (Delta,
Gamma, etc.) come immediately for free as well as various derivatives on
$\sigma$.

A weak point of the proposed method for the SABR model (but not for the Heston
model as this is shown in (Carr et al., 2021)) is taking into account
correlation between two Brownian motions of the model. We developed our
approach assuming zero correlation, and then showed how small correlations
could be taken into account. However, for strong correlations the asymptotic
approach is not feasible. Therefore, at the moment this remains to be an open
question how the GIT method could be applied to the $\lambda$-SABR model with
strong correlation. One idea is to switch dimensions, i.e. to construct the
GIT in the $\sigma$ space first and then it can be explicitly inverted as the
model is log-normal in $\sigma$, (Itkin et al., 2021). Then instead of the
LMVF equation in $t-\sigma$ variables we obtain a similar LMVF equation in
$t-x$ variables which now includes correlation. Solving it we can obtain the
option price which in this case is represented not via series, but is just
proportional to the solution of the LMVF equation. We will investigate this
approach in our future work.

## Acknowledgments

We are grateful to Peter Carr, Elisabeth Larsson, Alexander Lipton and
Fazlollah Soleymani for various discussions. Dmitry Muravey acknowledges
support by the Russian Science Foundation under the Grant number 20-68-47030.

## References

  * Abramowitz and Stegun (1964) M. Abramowitz and I. Stegun.  _Handbook of Mathematical Functions_.  Dover Publications, Inc., 1964.
  * Antonelli and Scarlatti (2009) F. Antonelli and S. Scarlatti.  Pricing options under stochastic volatility:a power series approach.  _Finance Stoch_ , 13:269–303, 2009.
  * Antonov et al. (2019) A. Antonov, M. Konikov, and M. Spector.  _Modern SABR Analytics_.  SpringerBriefs in Quantitative Finance. Springer, 2019.  ISBN 978-3-030-10656-0.
  * Assari et al. (2019) P. Assari, F. Asadi-Mehregan, and M. Dehghan.  On the numerical solution of Fredholm integral equations utilizing the local radial basis function method.  _International Journal of Computer Mathematics_ , 96(7):1416–1443, 2019.
  * Barger and Lorig (2017) W. Barger and M. Lorig.  Approximate pricing of european and barrier claims in a local-stochastic volatility setting.  _International Journal of Financial Engineering_ , 4(02n03):1750017, 2017.
  * Bateman and Erdélyi (1953) H. Bateman and A. Erdélyi.  _Higher Transcendental Functions_ , volume 1 of _Bateman Manuscript Project California Institute of Technology_.  McGraw-Hill, 1953.
  * Carr and Itkin (2021) P. Carr and A. Itkin.  Semi-closed form solutions for barrier and American options written on a time-dependent Ornstein-Uhlenbeck process.  _Journal of Derivatives_ , 29(1):9–26, 2021.
  * Carr and Sun (2007) P. Carr and J. Sun.  A new approach for option pricing under stochastic volatility.  _Review of Derivatives Research_ , 10:87–250, 2007.
  * Carr et al. (2020) P. Carr, A. Itkin, and D. Muravey.  Semi-closed form prices of barrier options in the time-dependent CEV and CIR models.  _Journal of Derivatives_ , 28(1):26–50, 2020\.
  * Carr et al. (2021) P. Carr, A. Itkin, and D. Muravey.  Semi-analytical pricing of barrier options in the time-dependent Heston model.  In preparation, 2021.
  * Fasshauer and McCourt (2012) G.E. Fasshauer and M.J. McCourt.  Stable evaluation of Gaussian Radial Basis function interpolants.  _SIAM Journal on Scientific Computing_ , 34:A737–A762, 2012\.
  * Fornberg et al. (2011) B. Fornberg, E. Larsson, and N. Flyer.  Stable computations with Gaussian Radial Basis functions.  _SIAM J. Sci. Comput._ , 33(2):869–892, 2011\.
  * Gorovoi and Linetsky (2004) V. Gorovoi and V Linetsky.  Black’s model of interest rates as options, eigenfunction expansions and japanese interest rates.  _Mathematical Finance_ , 14(1):49–78, 2004.
  * Gray and Pinsky (1992) A. Gray and M.A. Pinsky.  Computer graphics and a new Gibbs phenomenon for Fourier-Bessel series.  _Experimental Mathematics_ , 1(4):313–316, 1992\.
  * Hagan et al. (2002) P. Hagan, D. Kumar, A. Lesniewski A, and D Woodward.  Managing smile risk.  _Wilmott magazine_ , pages 84–108, September 2002.
  * Hagan et al. (2014) P.S. Hagan, D. Kumar, A. Lesniewski, and D. Woodward.  Arbitrage-free SABR.  _Wilmott_ , pages 60–75, January 2014.
  * Hagan et al. (2020) P.S Hagan, A. Lesniewski, and D.E Woodward.  Implied volatilities for mean reverting SABR models.  _Wilmott_ , pages 62–77, July 2020.
  * Hardy (1991) G.H. Hardy.  _Divergent Series_.  AMS Chelsea Publishing, 2 edition, 1991.  ISBN 978-0821826492.
  * Henry-Labordere (2005) P. Henry-Labordere.  A general asymptotic implied volatility for stochastic volatility models, 2005.  URL https://arxiv.org/pdf/cond-mat/0504317.pdf.
  * Itkin (2015) A. Itkin.  High-Order Splitting Methods for Forward PDEs and PIDEs.  _International Journal of Theoretical and Applied Finance_ , 18(5):1550031–1 —1550031–24, 2015.
  * Itkin (2017) A. Itkin.  _Pricing derivatives under Lévy models_.  Number 12 in Pseudo-Differential Operators. Birkhauser, Basel, 1 edition, 2017.
  * Itkin and Muravey (2021) A. Itkin and D. Muravey.  Semi-analytic pricing of double barrier options with time-dependent barriers and rebates at hit.  _Frontiers of Mathematical Finance_ , 1:1–36, 2021.
  * Itkin et al. (2020) A Itkin, A. Lipton, and D. Muravey.  From the Black-Karasinski to the Verhulst model to accommodate the unconventional Fed’s policy, June 2020.  URL https://arxiv.org/abs/2006.11976.
  * Itkin et al. (2021) A. Itkin, A. Lipton, and D. Muravey.  _Generalized Integral Transforms in Mathematical Finance_.  WSPC, Singapore, 2021.  ISBN 978-981-123-173-5.
  * Kartashov (2001) E.M. Kartashov.  _Analytical Methods in the Theory of Heat Conduction in Solids_.  Vysshaya Shkola, Moscow, 2001.
  * Kato et al. (2013) T. Kato, A. Takahashi, and T. Yamada.  An asymptotic expansion formula for Up-and-Out barrier option price under stochastic volatility model.  _JSIAM Letters_ , 5:17–20, 2013.
  * Larsson and Fornberg (2005) E. Larsson and B. Fornberg.  Theoretical and computational aspects of multivariate interpolation with increasingly flat radial basis functions.  _Comput. Math. Appl._ , 49(1):103–130, 2005.
  * Larsson et al. (2013) E. Larsson, E. Lehto, A. Heryudono, and B. Fornberg.  Stable computation of differentiation matrices and scattered node stencils based on Gaussian radial basis functions.  _SIAM J. Sci. Comput._ , 35(4):2096–2119, 2013\.
  * Lipton (2001) A. Lipton.  _Mathematical Methods For Foreign Exchange: A Financial Engineer’s Approach_.  World Scientific, 2001.
  * Mumford et al. (1983) D. Mumford, C. Musiliand M. Nori, E. Previato, and M. Stillman.  _Tata Lectures on Theta_.  Progress in Mathematics. Birkhäuser Boston, 1983.  ISBN 9780817631093.
  * Paige and Saunders (1975) C. Paige and M. Saunders.  Solution of sparse indefinite systems of linear equations.  _SIAM J. Numer. Anal._ , 12:617–629, 1975.
  * Platen (1997) E. Platen.  A non-linear stochastic volatility model.  Financial Mathematics Research Report FMRR005-97, Canberra: Center for Financial Mathematics, Australian National University, 1997.
  * Polyanin (2002) A.D. Polyanin.  _Handbook of linear partial differential equations for engineers and scientists_.  Chapman & Hall/CRC, 2002.
  * Revuz and Yor (1999) D. Revuz and M. Yor.  _Continuous Martingales and Brownian Motion_.  Springer, Berlin, Germany, 3rd edition, 1999.
  * Shreve (1992) S. Shreve.  Martingales and the theory of capital-asset pricing.  _Lecture Notes in Control and Information SCIENCES_ , 180:809–823, 1992.
  * Thakoor et al. (2019) N. Thakoor, D.Y. Tangman, and M.A. Bhuruth.  A spectral approach to pricing of arbitrage-free SABR discrete barrier options.  _Computational Economics_ , 54:1085–1111, 2019.
  * Van der Stoep et al. (2015) A.W. Van der Stoep, L.A. Grzelak, and C.W. Oosterlee.  The time-dependent FX-SABR model: Efficient calibration based on effective parameters.  _International Journal of Theoretical and Applied Finance_ , 18(6):1550042, 2015.
  * Watson (1966) G.N. Watson.  _A Treatise on the Theory of Bessel Functions_.  Cambridge University Press, Cambridge, UK, 2nd edition, 1966.
  * Yang et al. (2017) N. Yang, Y. Liu, and Z. Cui.  Pricing continuously monitored barrier options under the SABR model: A closed form approximation.  _JMSE_ , 2(2):116–131, 2017.
  * Zhang et al. (2014) H. Zhang, Y. Chen, and X. Nie.  Solving the linear integral equations based on radial basis function interpolation.  _Journal of Applied Mathematics_ , 2014:793582, 2014.