#  On a Numerical Solution of Structural model for the Probability of Default
under a Regime-Switching Synchronous-Jump Tempered Stable Lévy Model with
Desingularized Meshfree Collocation method

Davood Damircheli  Seyed-Mohammad-Mahdi Kazemi  Ali Foroush Bastani
Department and Organization, Mississippi State University,
P.O. Box 39759, Starkville, USA  Department of Financial Mathematics, Faculty
of Financial Sciences, Kharazmi University,
P.O. Box 15936-56311, Tehran, Iran  Department of Mathematics, Institute for
Advanced Studies in Basic Sciences,
P.O. Box 45195-1159, Zanjan, Iran

###### Abstract

In the paper [Hainaut, D. and Colwell, D.B., A structural model for credit
risk with switching processes and synchronous jumps, The European Journal of
Finance 22(11) (2016): 1040-1062], the authors exploit a synchronous-jump
regime-switching model to compute the default probability of a publicly traded
company. Here, we first generalize the proposed Lévy model to more general
setting of tempered stable processes recently introduced into the finance
literature. Based on the singularity of the resulting partial integro-
differential operator, we propose a general framework based on strictly
positive-definite functions to de-singularize the operator. We then analyze an
efficient meshfree collocation method based on radial basis functions to
approximate the solution of the corresponding system of partial integro-
differential equations arising from the structural credit risk model. We show
that under some regularity assumptions, our proposed method naturally de-
sinularizes the problem in the tempered stable case. Numerical results of
applying the method on some standard examples from the literature confirms the
accuracy of our theoretical results and numerical algorithm.

###### keywords:

Probability of Default; Structural Credit Risk Model; Strictly Positive-
Definite Functions; Radial Basis Function Collocation; Regime Switching;
Tempered Stable Lévy Process.

##  1 Introduction

In recent years, credit risk models have become an indispensable tool for
regulators to assess the performance of financial institutions, specially
banking and credit providing entities. Motivated by the proposals in Basel
accords111https://www.bis.org/bcbs/publ/d424.htm, accurate credit risk models
are now essential tools in predicting the default risk of loan portfolios,
evaluating the vulnerability of lender institutions and estimating reliable
economic capital levels for banks to remain solvent at a given confidence
level and time horizon [1]. They also provide useful tools for market
participants to identify, measure, monitor and control their exposure to
credit risk [2].

The two primary approaches in the literature to model the default risk of debt
obligations and to price credit risky securities are structural and reduced-
form (a.k.a intensity or statistical) models [3]. While the “reduced-form”
approach, abstracts away from the economic notion of bankruptcy and treats the
default event as an exogenous occurrence governed by a specific jump-diffusion
process, the structural models use the evolution of firms’ structural
variables, such as asset and debt values, to determine the time of default.
These last models are based on capital structure theory of Modigliani and
Miller [4] and option pricing theory of Black and Scholes [5] and Merton [6]
and rely heavily on diffusion processes to model the evolution of the firm’s
asset value process [7].

There have been some efforts in the literature to extend the Merton’s
structural model to more complex dynamical processes such as jump-diffusion
[8] or regime-switching diffusion processes [9]. The family of Markov-
modulated regime switching processes based on continuous-time Markov chains
has been extensively used and analyzed by academics and market practitioners
to model the variable economic conditions observed frequently in finance and
economics (see e.g. [10, 11] and the references therein).

In the case that we are confronted with a sudden synchronous jump in the asset
value process alongside the regime shifts, we obtain a process recently
introduced into the field of asset pricing (see Chourdakis [12] for a complete
account) with some option pricing applications (see e.g. [13, 14, 15]). It is
shown empirically that this regime-switching model and its extensions could
successfully be calibrated to a wide range of asymmetric volatility profiles
(see e.g. [16]).

In a recent contribution, Hainaut and Colwell [17] studied a synchronous jump
regime-switching version of Merton’s structural credit risk model to estimate
the probability of default and price corporate bonds. Based on the fact that
there is no closed-form solution for default probabilities and bond prices in
this setting, they developed a numerical approximation scheme based on the
Fourier space time-stepping method (see e.g. [18, 19]) which is only
applicable for Lévy processes with known characteristic functions.

In this paper, we extend the family of Lévy processes utilized in Hainaut and
Colwell [17] by considering the family of tempered stable distributions which
capture the behaviour of firm’s asset value process more accurately. Tempered
stable distributions and their associated processes are a class of models that
have attracted the attention of many researchers from applied probability and
stochastic analysis to physics and financial mathematics. They were introduced
in [20], where the associated Lévy processes are called “truncated Lévy
flights” and have been generalized by several authors [21]. Tempered stable
distributions form a six parameter family of infinitely divisible
distributions, which cover several well-known sub-classes such as Variance
Gamma (VG) distributions [22, 23], bilateral Gamma distributions [21] and CGMY
distributions [7]. Properties of tempered stable distributions have been
investigated, e.g., in [24, 25, 26, 27]. For financial modelling, they have
been applied, e.g. in [9, 28, 29, 3], see also the recent book [30].

Based on the fact that the system of partial integro-differential equations
(PIDEs) arising from these models contains singular integral terms arising
from the structure of their probability density functions are difficult to
treat numerically and these singularities adversely affect the accuracy and
convergence rate of the standard numerical methods (see e.g. [13, 31, 25, 32,
33, 34, 35, 36]), we need to employ special de-singularization techniques to
overcome these difficulties with minimum extra computational cost. Some
research studies have been done to address this issue in the literature. Cont
and his coworkers [37, 33, 29] have approximated the above-mentioned PIDEs
with a new singularly perturbed PIDE($\epsilon$) where small jumps near the
origin are estimated by an appropriate Brownian motion process. Exploiting the
Fourier transform technique is another alternative employed for PIDEs with and
without singularities [38]. However, another approach which is inspired form
physical intuition is to convert the original PIDEs into semi-parabolic PDEs
by eliminating the integral terms and applying a modern finite difference
method to solve the derived problem [35]. Last but not least, using a stable
method to approximate the solution of the PIDEs such that the numerical method
inherently handles the singularity is similar to the method proposed in this
paper.

In recent years, considerable attention has been centered around the
development of new numerical schemes under the heading of mesh-free or
meshless methods. They constitute powerful tools in the field of numerical
analysis with a wide range of applications from image processing to machine
learning [39]. Among these meshless methods, those which are based on
expansion in terms of radial basis functions (RBFs) provide a simple and
efficient framework to implement the meshfree idea and offer promising results
in real world applications. As a method to solve partial (integro-)
differential equations arising from a wide variety of applications, RBF-based
methods are used both in Galerkin and collocation modes with good accuracy and
stability properties (see e.g. [40, 41] and references therein).

In this paper, we propose and analyze an RBF-based collocation scheme based on
global multiquadratic basis functions with inherent capability to
desingularize the integral terms in a straightforward manner. Our contribution
could be considered as an extension of the method presented in Brummelhuis and
Chan [34] to the more general setting of tempered stable Lévy processes where
in addition, we have provided a firm theoretical basis on the working of the
proposed method. In order to validate the accuracy and efficiency of this
method, we consider the default probability estimation for three different
firms previously studied in [17] now based on a tempered stable Lévy
distribution as the firm’s asset value process. Our findings confirm the
capabilities of the new method as a reliable and accurate scheme to handle the
problem complexities.

The outline of the article is as follows. In Section 2, we introduce the basic
ingredients of our proposed model for the asset value process including a
brief account of the hidden Markov processes underlying the regime change
behaviour of the market and encoding the information about changing economic
conditions. We then provide some general comments about structural models to
compute the default probability of a zero coupon bond. Finally, we describe
the details of the underlying Markov-modulated switching Lévy processes with
synchronous jumps to model the asset value process of the firm under the risk-
neutral measure. In Section 3, we review the structural model for corporate
debts and the features of the firm’s asset value process modeled as an
exponential switching Lévy process. Among the numerous Lévy processes Next,
the system of partial integro-differential equations (PIDEs) driving the
probabilities of default is built. Sections 4 present the numerical method
based on the proposed radial basis function collocation. In Section 5, some
numerical experiments are conducted to illustrate the efficiency and accuracy
of the proposed method. We conclude the paper by pointing out some future
research directions.

##  2 Proposed Model for the Value of Firm’s Assets

In order to model the firm’s asset value process, we employ a Markov-modulated
exponential Lévy process having synchronous jumps with regime shifts. Before
presenting the details of this process, we first remind the reader of some
general properties of continuous-time Markov processes. In this respect, we
assume that the economic state at time $t$ is modeled by a continuous-time
hidden Markov process, $\alpha_{t}$, with values from the set
$\mathcal{H}=\\{1,2,\cdots,H\\}$, each representing an economic state or
regime. Let $Q=(q_{ij})_{H\times H}$ denote the generator matrix of
$\alpha_{t}$ where the off-diagonal entries are non-negative and
$\sum_{j=1}^{H}q_{ij}=0,$ for all $i\in\mathcal{H}$ (see e.g. [42]).
Discretizing the time variable with a sufficiently small step-size, $\Delta
t$, we could show that $q_{i,j}\Delta t$ is the probability to switch from
state $i$ to state $j$, for $i\neq j$. We also define the transition
probability matrix of the process, denoted by $P(t,s)$ in the form

$P(t,s)={\rm e}^{Q(s-t)},\quad s\geq t.$ (2.1)

We denote the elements of $P(t,s)$ by $p_{ij}(t,s)$ and interpret them as
being probabilities of switching from state $i$ at time $t$, to state $j$ at
time $s$.

###  2.1 Switching Lévy Processes with Synchronous Jumps

In the remainder, we consider an exponential Lévy model for the asset value
process of the form

$V_{t}=V_{0}e^{X_{t}},$

in which $X_{t}$ is a Markov-modulated Lévy process222A Lévy process is a
stochastic process with stationary and independent increments which is
continuous in probability (see e.g. Papapantoleon [43]). with $X_{0}=0$. More
precisely, $X_{t}$ is a piecewise Lévy process depending upon the state,
$\alpha_{t}$ of an economy with simultaneous jumps at regime switching times.
Dynamics of the process $X_{t}$ could be written as

$dX_{t}=dX^{\alpha_{t}}_{t}+J_{\alpha_{t}^{-},\alpha_{t}},$ (2.2)

where $X_{t}$ has jumps represented by $J_{\alpha_{t}^{-},\alpha_{t}}$ when
regime switches at time $t$. The processes $X^{j}_{t}$, $j=1,\cdots H$ are
independent Lévy processes with a Lévy-Itô decomposition of the form

$\mathrm{d}X_{t}^{j}=\mu_{j}\mathrm{d}t+\sigma_{j}\mathrm{d}W_{t}^{j}+\int_{|z|>1}zJ_{X^{j}}(\mathrm{d}t,\mathrm{d}z)+\int_{|z|\leq
1}z\left(J_{X^{j}}(\mathrm{d}t,\mathrm{d}z)-\nu(j,\mathrm{d}z)\mathrm{d}t\right),$
(2.3)

in which $W^{j}_{t}$ is a standard Brownian motion on the underlying
probability space and $J_{X^{j}}(t,z)$ is a jump process of intensity
$\nu(j,\cdot)$ (different from $J_{\alpha_{t}^{-},\alpha_{t}}$), the Lévy
measure of $X_{t}$ in state $j$. The triplet $(\mu_{j},\sigma_{j},\nu(j,z))$
uniquely determines the characteristic function of $X^{j}_{t}$ by the Lévy-
Khintchine theorem (see e.g. [29, 20]) stating that

${\mathbb{E}}\left({\mathrm{e}}^{{\mathrm{i}}\omega
X_{t}^{j}}\right)={\mathrm{e}}^{t\Psi_{j}(\omega)},$ (2.4)

in which the characteristic exponent of the Lévy process is given by

$\Psi_{j}(w)={\mathrm{i}}\mu_{j}\omega-\frac{\omega^{2}\sigma_{j}^{2}}{2}+\int_{\mathbb{R}}\left({\mathrm{e}}^{{\omega}z}-1-{\mathrm{i}}\omega
z\mathbbm{1}_{|z|<1}\right)\nu(j,dz).$ (2.5)

In the above expressions, $\mu_{j}$ is the drift, $\sigma_{j}>0$ is the
diffusion and $\nu(j,\cdot)$ is the Lévy measure on $\mathbb{R}$ satisfying
$\nu(j,\\{0\\})=0$ and $\int_{\mathbb{R}}(1\wedge|z|^{2})\nu(j,dz)<\infty$
(see e.g. [20] for more details on Lévy processes).

The finite-activity Lévy density of synchronous jumps $J_{i,j}$ is given by
(see [17])

$\mu(i,j,z)=\left\\{{\begin{array}[]{*{20}{c}}{{\eta_{i,j}}{{\rm{e}}^{-{\eta_{i,j}}\left|z\right|}},}&{z\in{\mathbb{R}^{{\rm{sign}}({J_{i,j}})}},}\\\
{0,}&{{\rm{otherwise,}}}\end{array}}\right.$ (2.6)

which can take positive, negative or real values with a characteristic
function of the form

$\displaystyle{\Theta^{ij}}(u)$
$\displaystyle=\mathbb{E}\left(\mathrm{e}^{\mathrm{i}uJ_{i,j}}\right)=\int_{-\infty}^{+\infty}\mathrm{e}^{\mathrm{i}uz}\mu(i,j,z)\mathrm{d}z=\frac{{{\eta_{i,j}}}}{{{\eta_{i,j}}-{\rm{i}}u{\rm{sign}}({J_{i,j}})}}.$
(2.7)

In the next subsection, we introduce the specific Lévy process used in this
paper along with its Lévy measure.

###  2.2 Generalized Tempered Stable Lévy Processes

Consider a Lévy process, $X^{j}_{t}$, specified by it’s characteristic
triplet, $(\mu_{j},\sigma_{j},\nu(j,z))$. The process is called of finite
activity when $\nu(j,\mathbb{R})=\int_{\mathbb{R}}\nu(j,dz)<\infty$ and of
infinite activity otherwise (see e.g. [29]). We know also that a Lévy process
can be decomposed in the form
$X_{t}^{j}=\mu_{j}t+\sigma_{j}W^{j}_{t}+J^{j}_{X}$ where $J^{j}_{X}$ describes
a jump process that its jumps can be finite or infinite.

Let the family of functions, $\bf{W}$, to contain those
$w:(0,\infty)\rightarrow[0,\infty)$ which are continuous and also satisfy the
following conditions:

  * •

They are decreasing,

  * •

$w(0^{+})=1$,

  * •

$\mathop{\lim}\limits_{x\to\infty}{z^{n}}w(z)=0$, $\forall n\in\mathbb{N}$.

######  Definition 2.1.

A Lévy process is of generalized tempered-stable (GTS) type if its Lévy
measure, $\nu(j,\cdot)$, could be expressed as

${\nu}(j,dz):=\left({\frac{{{C^{j}_{+}}w({\beta^{j}_{+}}z)}}{{{z^{{1}+{\alpha^{j}_{+}}}}}}{\mathbbm{1}_{z>0}}(z)+\frac{{{C^{j}_{-}}w({\beta^{j}_{-}}\left|z\right|)}}{{{{\left|z\right|}^{1+{\alpha^{j}_{-}}}}}}{\mathbbm{1}_{z<0}}(z)}\right)dz,$
(2.8)

in which $w\in\bf{W}$ and the parameters $\alpha^{j}_{\pm}<2$ and
$C^{j}_{\pm},\beta^{j}_{\pm}>0$ are chosen such that
$|z|^{n}{C_{\pm}^{j}}w({\beta^{j}_{\pm}}z)$ is bounded on $\mathbb{R}$ for
each $n\in{\mathbb{N}}$.

######  Remark 2.2.

It could be shown that for $\alpha^{j}_{\pm}<2$, we have

$\int_{|z|>1}|z|^{n}\nu(j,dz)<\infty,\quad\forall\ n\geq 0.$ (2.9)

Indeed, the relation (2.9) is a necessary and sufficient condition for moments
of a GTS Lévy process to exist [29].

Table 1: Classification of Lévy processes based on their Lévy–Khintchin
representation. In this table, we introduce the parameters
$\tilde{\mu_{j}}=\mu_{j}-\int_{|z|<1}z\nu(j,dz)$ and
$\hat{\mu}_{j}=\tilde{\mu}_{j}-\int_{\mathbb{R}}\nu(j,dz)$ (For details see
[20, 29, 21, 34]).

.  Different types of Lévy jump-diffusion processes and their Lévy–Khintchin
representation Activity of $X^{j}_{t}$ Variation of $J^{j}_{X}$ Characteristic
exponent $\Psi_{j}(w)$ Parameters for GTStnote:robots-r1 Financial Models
Financial Refs. Finite Finitetnote:robots-r2
${\mathrm{i}}\hat{\mu}_{j}w+\dfrac{\sigma_{j}^{2}w^{2}}{2}+\int_{\mathbb{R}}e^{{\mathrm{i}}\omega
z}\nu(j,dz)$ $\alpha^{j}_{\pm}<0$ Merton and Kou [44] and [45] Infinite Finite
${\mathrm{i}}{\tilde{\mu}_{j}}w+\dfrac{{\sigma_{j}}^{2}w^{2}}{2}+\int_{\mathbb{R}}\left(e^{{\mathrm{i}}\omega
z}-1\right)\nu(j,dz)$ $\alpha^{j}_{\pm}\in[0,1)$ VG and CGMY [23] and [22]
Infinite
${\mathrm{i}}\mu_{j}w+\dfrac{\sigma^{2}w^{2}}{2}+\int_{\mathbb{R}}\left(e^{{\mathrm{i}}\omega
z}-1-{\mathrm{i}}wz\mathbbm{1}_{|z|\leq 1}\right)\nu(j,dz)$
$\alpha^{j}_{\pm}\in[1,2)$ CGMY [22] a Note that $M_{GTS}$ is not a Lévy
density when $\alpha^{j}_{+}\geq 2$ or $\alpha^{j}_{-}\geq 2$. b If a Lévy
jump-diffusion process is of finite activity, then it has also jumps of finite
variation.

######  Remark 2.3.

The following particular cases are known in the literature:

  * •

$C^{j}_{+}=C^{j}_{-}$ and $\alpha^{j}_{+}=\alpha^{j}_{-}=0$ is a Variance
Gamma (VG) distribution, see [24, 23];

  * •

$C^{j}_{+}=C^{j}_{-}$ and $\alpha^{j}_{+}=\alpha^{j}_{-}$ is a CGMY-
distribution, see [22], also called classical tempered stable distribution in
[28];

  * •

$\alpha^{j}_{+}=\alpha^{j}_{-}$ is a KoBol distribution, see [46];

  * •

$\alpha^{j}_{+}=\alpha^{j}_{-}$ and $\beta^{j}_{+}=\beta^{j}_{-}$ is the
infinitely divisible distribution associated to a truncated Lévy flight, see
[47];

  * •

$\alpha^{j}_{+}=\alpha^{j}_{-}=0$ is a bilateral Gamma distribution, see [30].

###  2.3 The Structural Model of Default Risk

Let us consider a firm whose capital structure consists of debt and equity
where the debt component is issued as a single zero-coupon bond, $B$ with
principal value $L$. We also assume that the firm can only default at the
maturity of debt, $T$. Assume also a risk-neutral probability space
$(\Omega,\mathcal{F},\mathbb{Q})$ equipped with some filtration,
$\mathcal{F}_{t}$. On this filtration is defined the firm’s value process,
$V_{t}$, that represents the total value of the firm’s assets. The probability
measure $\mathbb{Q}$ is taken to be a risk neutral measure but the following
results are also valid under the real historical measure, $\mathbb{P}$.

We define the default event to occur only if the total value of the firm’s
assets at the maturity is not sufficient to cover the loan’s redemption:

$\textmd{Default Event}=\\{V_{T}<L\\},$ (2.10)

and the probability of default is denoted by

$D(t,T)=\mathbb{P}(V_{T}\leq L|{\cal F}_{t}).$ (2.11)

In case of bankruptcy, a fraction $R(\alpha_{T})$, called the recovery rate,
of asset value is repaid to debt-holders.

In this work, as in Guo et al. [48], the debt recovers at a different
magnitude depending here upon the economic regime, $\alpha(t)$. If $r$ is the
constant risk-free rate of interest, then the time $t$ price of a defaultable
bond issued by the firm will be obtained as

$B(t,T)={\mathbb{E}}({\rm
e}^{-r(T-t)}(R(\alpha_{T})V_{T}\mathbbm{1}_{V_{T}<L}+L\mathbbm{1}_{VT\geq
L})|{\mathcal{F}}_{t}).$ (2.12)

We will present the approach used to compute the probabilities of default and
the defaultable bond prices in the reminder.

##  3 Derivation of the PIDEs for Default Probabilities

In order to derive the PIDEs corresponding to the default probability of a
given firm, there exists two main approaches which are usually exploited and
give equal results:

  1. 1.

Using Itô’s formula for general semi-martingales (or infinitesimal generator)
[29];

  2. 2.

Using the Lévy–Khintchine formula and the tempered Fourier transform (see
Proposition (1.9) of [26]).

The second approach is more useful in characterizing some intrinsic properties
of the integro-differential operator of the problem which will be needed in
later parts of the paper.

The probability of default of a firm in state $\alpha_{t}$ and time $t$ under
the measure $\mathbb{P}$ is defined as

$\displaystyle D(t,T,X_{t},{\alpha_{t}})$
$\displaystyle=\mathbb{P}\left(V_{T}\leq
L|\mathcal{F}_{t}\right)=\mathbb{E}\left(\mathbbm{1}_{V_{T}\leq
L}|\mathcal{F}_{t}\right),$ (3.1)

in which $D\left(t,T,x,j\right)$ is the probability of default at regime $j$
when the state process $X_{t}$ takes the value $x$. It could be shown (see
e.g. Hainaut and Colwell [11]) that $D\left(t,T,x,j\right)$ for
$j=1,2,\cdots,H$ satisfy the following system of partial integro-differential
equations

$\left\\{\begin{array}[]{ll}\dfrac{\partial}{\partial
t}D(t,T,x,j)+\mathcal{L}D(t,T,x,j)=0,\hskip 56.9055ptj=1\ldots H,\\\ \\\
D(t,T,x,j)=\mathbbm{1}_{x<\ln\left(L/V_{0}\right)},\hskip 95.3169ptj=1\ldots
H,\\\ \\\ \lim_{x\rightarrow-\infty}D(t,T,x,j)=1,\hskip 93.89418ptj=1\ldots
H,\\\ \\\ \lim_{x\rightarrow+\infty}D(t,T,x,j)=0,\hskip 93.89418ptj=1\ldots
H,\end{array}\right.$ (3.2)

in which $\mathcal{L}$ is the infinitesimal generator of default probabilities
of the form (see e.g. [17, 29])

$\displaystyle\mathcal{L}D(t,T,x,j)=$
$\displaystyle\mu_{j}\frac{\partial}{\partial
x}D(t,T,x,j)+\frac{\sigma_{j}^{2}}{2}\frac{\partial}{\partial
x^{2}}D(t,T,x,j)$ (3.3) $\displaystyle+\sum_{k\neq
j}q_{j,k}\int_{\mathbb{R}\backslash\\{0\\}}D(t,T,x+z,k)-D(t,T,x,j)\mu(j,k,z)\mathrm{d}z$
$\displaystyle+\int_{\mathbb{R}\backslash\\{0\\}}\left(D(t,T,x+z,j)-D(t,T,x,j)-z\mathbbm{1}_{|z|\leq
1}\frac{\partial}{\partial x}D(t,T,x,j)\right)v(j,\mathrm{d}z).$

Due to the relation $\sum_{i=1}^{H}q_{j,i}=0$, the generator can be written as

$\displaystyle\mathcal{L}D(t,T,x,j)=$
$\displaystyle\mu_{j}\frac{\partial}{\partial
x}D(t,T,x,j)+\frac{\sigma_{j}^{2}}{2}\frac{\partial}{\partial
x^{2}}D(t,T,x,j)$ (3.4) $\displaystyle+q_{j,j}D(t,T,x,j)+\sum_{k\neq
j}q_{j,k}\int_{\mathbb{R}\backslash\\{0\\}}D(t,T,x+z,k)\mu(j,k,z)\mathrm{d}z$
$\displaystyle+\int_{\mathbb{R}\backslash\\{0\\}}\left(D(t,T,x+z,j)-D(t,T,x,j)-z\mathbbm{1}_{|z|\leq
1}\frac{\partial}{\partial x}D(t,T,x,j)\right)v(j,\mathrm{d}z).$

In order to make the presentation and implementation of the numerical meshfree
scheme in the next section simple, we transform the PIDEs (3.2) forward in
time. To this end and using a common change of variables of the form
$t=T-\tau$ and thus $u(\tau,x,j)=D(t,T,x,j)$ where $T$ is removed from $u$ for
simplicity of notation, the problem (3.2) will be transformed into

$\left\\{\begin{array}[]{ll}\dfrac{\partial}{\partial\tau}u(\tau,x,j)=\mathcal{L}u(\tau,x,j),\hskip
62.59596ptj=1\ldots H,\\\ \\\
u(0,x,j)=\mathbbm{1}_{x<\ln\left(L/V_{0}\right)},\hskip 73.97733ptj=1\ldots
H,\\\ \\\ \lim_{x\rightarrow-\infty}u(\tau,x,j)=1,\hskip 71.13188ptj=1\ldots
H,\\\ \\\ \lim_{x\rightarrow+\infty}u(\tau,x,j)=0,\hskip 71.13188ptj=1\ldots
H.\end{array}\right.$ (3.5)

Developing a stable numerical scheme to solve Eq. (3.2) is a necessity as the
analytical solution of the PIDE system is not available in closed form. In the
sequel, we will propose our approach to approximate the solution of these
systems of partial integro-differential equations.

##  4 De-Singularized Meshfree Approximation

The important point to notice in confronting with the infinite activity
processes is the singularity of the Lévy measure $\nu(j,dz)$ near the origin
due to small jumps. This Lévy measure for different processes behaves
differently near zero. For instance, the Lévy measure for the CGMY process (or
equally the GTS process with parameters $C^{j}_{+}=C^{j}_{-}$ and
$\alpha^{j}_{\pm}>0$) approaches infinity much faster than the variance gamma
(VG) process [22]. For this reason we should have the highest possible order
of approximation for the integrand near zero.

Indeed, as we move in Table 1 from top rows to bottom, more complex Lévy
processes appear due to the limitations which appear in separating the three
integrands (see e.g. the third column in this Table). So we will need more
complex numerical algorithms which will treat the singular behaviour more
efficiently (see e.g. [13, 31, 25, 32, 33, 34, 35]).

Based on the above discussion, various numerical algorithms (mainly based on
finite differences) are proposed based on different kinds of integral terms
appearing in PIDEs. Some researchers in the field have tried to offer
comprehensive algorithms to tackle different Lévy densities both with or
without singularity. Following is the list of outstanding works in this
subject:

  * •

Approximating small jumps near the origin using an appropriate approximated
Brownian motion which will result in a new PIDE. Since this PIDE is a function
of $\epsilon$, PIDE($\epsilon$) is used to refer to this derived PIDE. When
$\epsilon$ goes to $0$, it is proved that the solution of PIDE($\epsilon$)
will converge to the solution of original PIDE ([37, 33, 29]). Applying this
approach to the system (3.2) will lead us to the following PIDE($\epsilon$)

$\frac{\partial}{\partial\tau}u+\left(\mu_{j}-\omega_{j}(\varepsilon)\right)\frac{\partial}{\partial
x}u+\left(\frac{1}{2}\left(\sigma_{j}^{2}+\sigma_{j}(\varepsilon)\right)\frac{\partial^{2}u}{\partial
x^{2}}+q_{j,j}u\right)\\\ +\sum_{k\neq
j}q_{j,k}\int_{{\mathbb{R}/\\{0\\}}}u(x+z)\mu(j,k,dz)+\int_{|z|>\varepsilon}(u(\tau,j,x+z)-u(\tau,j,x))\nu(j,dz)=0,$
(4.1)

where the factors $\omega_{j}(\varepsilon)$ and $\sigma_{j}(\varepsilon)$ have
been introduced in [29] (Chapter 12 ).

  * •

An alternative methodology to handle such singularities is the famous Fourier
Transform technique in which the PIDE (with or without singularity) is
transformed into a system of ODEs for which efficient numerical approximation
algorithms exist. This method is usually used to benchmark the financial
problems (for more details, we refer the reader to [38]).

  * •

Last but not least, we will mention the method used by Itkin et al. [35] in
which by eliminating the integral part form original PIDEs, we deal with a
semi-parabolic PDE. where modern finite difference methods are used to
approximate the solution of this semi-parabolic PDE. It is worth noticing that
these family of PDEs have a complex structure so that the standard numerical
algorithms are not able to approximate their solutions. For details refer to

Our proposed approach is to use a meshless collocation method based on RBFs of
infinite or finite smoothness. Similar to Brummelhuis and Chan (see Section
3.3 of [34]), we show that smoothness properties of radial basis functions
such as MQ can overcome the singularity of a wide range of important Lévy
densities.

###  4.1 Problem Discretization

In this section, we describe a general numerical scheme based on RBF
collocation to numerically solve the transformed linear partial integro-
differential equation formulation (3.2) which could be represented as

$U_{\tau}={\cal L}U,$ (4.2)

in which

$U(\tau,x):=[u(\tau,x,1),u(\tau,x,2),\cdots,u(\tau,x,H)]^{T}.$ (4.3)

Many numerical schemes for time-dependent PDEs separate the discretization of
time and space variables into distinct phases and develop the theory by
assuming one discretization (outer discretization) to be carried out first
leading to a so-called semi-discrete problem. After investigating the thus
arising type of problem, one continues to perform the second discretization
(inner discretization), ending up with a fully discrete scheme (see e.g.
[49]).

We follow this tradition by first discretizing the equation (3.5) in time and
then using Newton’s method for linearizing the time-discrete problem at the
PDE level. Such a technique transforms the nonlinear stationary PDE at each
time level into a sequence of linear PDEs which could now be solved using RBF
collocation based on multiquadric (MQ) radial basis functions (see [50] and
[51] for more details).

###  4.2 Semi-Discretization in Time

Given an equally spaced time grid $0=\tau_{0}<\tau_{1}<...<\tau_{N}=T$ which
subdivides the interval $[0,T]$ into $N$ sub-intervals of the form
$[\tau_{n},\tau_{n+1}],~{}n=0,1,2,\cdots,N-1$ with $\tau_{n}=n\Delta\tau$ and
$\Delta\tau=T/N$, we first employ the temporal semi-discretization
$U(\tau_{n},x)\approx U_{n}(x)$

$\frac{{U_{n+1}(x)-U_{n}(x)}}{{\Delta\tau}}=\theta{\cal
L}[{U}]\big{|}_{\tau=\tau_{n}}+(1-\theta){\cal
L}[{U}]\big{|}_{\tau=\tau_{n+1}},$ (4.4)

to arrive at

$[1-\Delta\tau(1-\theta){\cal L}]{U}_{n+1}(x)=[1+\Delta\tau\theta{\cal
L}]{U}_{n}(x),\quad n=0,1,\cdots,N-1,$ (4.5)

in which $\theta\in[0,1]$ is the implicitness level of the scheme.

This last expression could be re-phrased as a linear elliptic partial
differential equation of the form

$\displaystyle{\cal L}_{\Delta\tau}U_{n+1}(x)={F}_{n},\quad n=0,\cdots,N-1,$
(4.6)

with

${\cal L}_{\Delta\tau}:=1-\Delta\tau(1-\theta){\cal
L},\quad{F}_{n}:=[1+\Delta\tau\theta{\cal L}]U_{n}(x).$ (4.7)

###  4.3 Meshfree RBF Collocation

In order to give a rough idea of the radial basis function collocation, we
first need the following concepts:

######  Definition 4.4.

A multivariate function $\Phi:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ is
called radially symmetric or isotropic on ${\mathbb{R}}^{d}$, if $\Phi({\bf
x})=\Phi({\bf y})$ whenever $\|{\bf x}\|=\|{\bf y}\|$. A radially symmetric
function $\Phi(\cdot)$ could be represented as
$\Phi(\bf{x})=\phi(\|\bf{x}\|)$, $\bf{x}\in$ ${\mathbb{R}}^{d}$ for some
univariate function, $\phi:[0,\infty)\rightarrow\mathbb{R}$. Note that the
norms will be the usual Euclidean 2-norm, defined by $\|{\bf
x}\|_{2}:=\sqrt{{\bf x}^{T}{\bf x}}$.

Before presenting the method, we first truncate the infinite spatial domain of
the problem into a finite sub-domain of the form $I=[x_{min},x_{max}]$ and
then consider the projection of $U(\tau,x)$ onto a finite dimensional space as
a linear combination of radial basis function of the form

$U_{n}(x)=\sum_{j=1}^{N_{x}}\Upsilon_{j}^{(n)}\phi\left(\left\|x-x_{j}\right\|\right)\equiv\sum_{j=1}^{N_{x}}\Upsilon_{j}^{(n)}\Phi_{j}(x),$
(4.8)

in which, $x_{j}$ for $j=1,\cdots,N_{s}$ is a set of scattered data points in
the sub domain $I$ and
$\Upsilon_{j}^{(n)}:=[\upsilon_{j1}^{(n)},\upsilon_{j2}^{(n)},\cdots,\upsilon_{jH}^{(n)}]^{T}$.
Now, if we use $\theta=0$ and Substitute (4.8) into (4.6), we obtain

$\sum_{j=1}^{N_{x}}\Upsilon_{j}^{(n+1)}\Phi_{j}(x)-\Delta\tau\sum_{j=1}^{N_{x}}\Upsilon_{j}^{(n+1)}\mathcal{L}\Phi_{j}(x)=U_{n}(x).$
(4.9)

So we have

$(\Phi-\Delta\tau\Phi_{\mathcal{L}})\Upsilon^{(n+1)}={\bf F}^{(n)}.$ (4.10)

where

$\displaystyle\Phi=\left[\Phi_{j}(x_{i})\right]_{i,j=1,\cdots,N_{x}},$ (4.11)
$\displaystyle\Phi_{\mathcal{L}}=\left[\mathcal{L}\Phi_{j}(x_{i})\right]_{i,j=1,\cdots,N_{x}},$
(4.12) $\displaystyle{\bf
F}^{(n)}:=[U_{n}(x_{1}),U_{n}(x_{2}),\cdots,U_{n}(x_{N_{x}})],$ (4.13)
$\displaystyle\Upsilon^{(n+1)}:=\left({\Upsilon_{1}^{(n+1)}}^{T},{\Upsilon_{2}^{(n+1)}}^{T},\ldots,{\Upsilon_{N_{x}}^{(n+1)}}^{T}\right)^{T}.$
(4.14)

It is worth mentioning that mega matrix $\Phi_{\mathcal{L}}$ has a block
structural where off-diagonal blocks entangles data from different regimes
(see for more details [13]).

###  4.4 Some Basic Facts about Radial Basis Function Interpolation

Let us recall some fundamental properties of RBFs before presenting the core
theorem of this section.

###  4.5 De-Singularization Procedure

The first part of this subsection is to review some of the necessary concepts
and theorems which we will apply to prove the main theorem of this article.

######  Definition 4.5.

(Bochner [52]) A radially symmetric multivariate function
$\Phi:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}$ is called positive definite, if
for any finite $k\in{\mathbb{N}}$ pairwise different points ${\bf x}_{1},{\bf
x}_{2},\cdots,{\bf x}_{k}$ in ${\mathbb{R}}^{d}$ and for each ${\bf
c}=[c_{1},c_{2},\cdots,c_{k}]^{T}\in{\mathbb{R}}^{k}$, we have

$\sum_{i=1}^{k}\sum_{j=1}^{k}c_{i}c_{j}\Phi({\bf x}_{i}-{\bf x}_{j})\geq 0.$
(4.15)

The function $\Phi$ is called strictly positive definite on
${\mathbb{R}}^{d}$, if the quadratic form in the left hand side of (4.15) is
zero only for ${\bf c}={\bf 0}$. We also call a univariate function
$\phi:[0,\infty)\rightarrow\mathbb{R}$ strictly positive definite on
$\mathbb{R}^{d}$ (abbreviated as $\phi\in{\bf SPD}_{d}$), if the corresponding
radially symmetric multivariate function defined by
$\Phi(\bf{x}):=\phi(\|\bf{x}\|)$, $\bf{x}\in$ ${\mathbb{R}}^{d}$ is strictly
positive definite.

In order to prove our main result, we need the following properties of ${\bf
SPD}_{d}$ functions which play a key role in our development.

######  Theorem 4.6.

(Fasshauer [39]) Every ${\bf SPD}_{d}$ function, $\Phi(\cdot)$, satisfies the
following properties:

  * •

$\Phi(0)\geq 0$;

  * •

$|\Phi(x)|\leq\Phi(0)$;

  * •

$\Phi$ is a positive function.

######  Theorem 4.7.

(Buescu and Paixao [53] and Massa and et al. [54]) Let
$\Phi:\mathbb{R}^{m}\rightarrow\mathbb{R}$ be a positive definite function and
suppose $\Phi$ is of class $C^{2n}$ (resp. $C^{\infty}$) in some neighborhood
of the origin for some positive integer $n$. Then

  * •

$\Phi\in C^{2n}(\mathbb{R}^{m})$ (resp. $C^{\infty}(\mathbb{R}^{m}))$,

  * •

${\left|{{{\rm{D}}^{\alpha+\beta}}\Phi(x)}\right|^{2}}\leq{(-1)^{{m_{1}}+{m_{2}}}}{{\rm{D}}^{2\alpha}}\Phi(0){{\rm{D}}^{2\beta}}\Phi(0),\quad\forall
x\in\mathbb{R}^{m},\quad|\alpha|,|\beta|\leq n$,

  * •

In the special case where $m=1$, we have $\Phi\in C^{2n}(\mathbb{R})$ and for
all integers $m_{1}$, $m_{2}$ with $0\leq m_{i}\leq n$, $i=1,2$ and every
$x\in\mathbb{R}$ we have

${\left|{{\Phi^{({m_{1}}+{m_{2}})}}(x)}\right|^{2}}\leq{(-1)^{{m_{1}}+{m_{2}}}}{\Phi^{(2{m_{1}})}}(0){\Phi^{(2{m_{2}})}}(0).$

######  Remark 4.8.

Based on the above results, we have the following:

  * (i)

Every ${\bf SPD}_{d}$ function, $\Phi$, is bounded above (see the Theorem
4.6).

  * (ii)

Based on above assumptions, $\Phi\in{\bf SPD}_{d}\cap C^{2n}(\mathbb{R})$ for
$n\geq 2$ ($\Phi$ is $2n$-times differentiability on $\mathbb{R}$).

  * (iii)

Based on the Theorem 4.7, if we define
${M_{k}}:=\mathop{\max}\limits_{x\in\mathbb{R}}\\{|{{\Phi^{(k)}}(x)}|\\}$ for
$0\leq k\leq 2n$, then $M_{k}<\infty$.

######  Theorem 4.9.

Let $\Phi:\mathbb{R}\rightarrow\mathbb{R}$ be a real-valued strictly positive
definite function and assume that $\Phi$ is of class $C^{2n}$ in some
neighborhood of the origin for some positive integer $n\geq 2$. Then, we have
$\Phi\in C^{2n}(\mathbb{R})$. Also for any tempered stable Lévy measure
$\nu(i,\mathrm{d}z)$ with $\alpha^{j}_{\pm}<2$ defined in 2.8, the integral

$I_{v}(x):=\int_{\mathbb{R}}\left(\Phi(x+z)-\Phi(x)-z\mathbbm{1}_{|z|\leq
1}\frac{\partial}{\partial x}\Phi(x)\right)\nu(j,\mathrm{d}z),$ (4.16)

is finite for all $x\in\mathbb{R}$ and without singularity at the origin.

###### Proof.

For arbitrary constants $0<\varepsilon<1$ and $b\gg 1$, the improper integral
$I_{v}(\cdot)$ could be written as

$I_{v}=\mathop{\lim}\limits_{b\to\infty}\left(I^{(1)}_{v}+I^{(4)}_{v}\right)+\mathop{\lim}\limits_{\varepsilon\to
0^{+}}\left(I^{(2)}_{v}+I^{(3)}_{v}\right),$

where

$\displaystyle I^{(1)}_{v}$
$\displaystyle=\int_{-b}^{-1}\left(\Phi(x+z)-\Phi(x)\right)\frac{{{C^{j}_{-}}w({\beta^{j}_{-}}\left|z\right|)}}{{{{\left|z\right|}^{1+{\alpha^{j}_{-}}}}}}\mathrm{d}z,$
$\displaystyle I^{(2)}_{v}$ $\displaystyle=\int_{-1}^{-\varepsilon}\
\left(\Phi(x+z)-\Phi(x)-z\mathbbm{1}_{|z|\leq 1}\frac{\partial}{\partial
x}\Phi(x)\right)\frac{{{C^{j}_{-}}w({\beta^{j}_{-}}\left|z\right|)}}{{{{\left|z\right|}^{1+{\alpha^{j}_{-}}}}}}\mathrm{d}z,$
$\displaystyle I^{(3)}_{v}$ $\displaystyle=\int_{\varepsilon}^{1}\
\left(\Phi(x+z)-\Phi(x)-z\mathbbm{1}_{|z|\leq 1}\frac{\partial}{\partial
x}\Phi(x)\right)\frac{{{C^{j}_{+}}w({\beta^{j}_{+}}z)}}{{{z^{{1}+{\alpha^{j}_{+}}}}}}\mathrm{d}z,$
$\displaystyle I^{(4)}_{v}$
$\displaystyle=\int_{1}^{b}\left(\Phi(x+z)-\Phi(x)\right)\frac{{{C^{j}_{+}}w({\beta^{j}_{+}}z)}}{{{z^{{1}+{\alpha^{j}_{+}}}}}}\mathrm{d}z.$

According to the comparison test for improper integrals (see Chapter 12 in
[55]), it is easy to show that the improper integrals of the first kind (i.e.
$I_{v}^{(1)}$ and $I_{v}^{(4)}$) are absolutely convergent as $b$ tends to
infinity for arbitrary $\alpha^{j}_{\pm}<2$, because we can obtain the
following upper bounds for the integrals $I^{(1)}_{v}$ and $I^{(4)}_{v}$ using
Theorem 4.6 and Definition 2.1:

$\displaystyle\left|I^{(1)}_{v}\right|$
$\displaystyle\leq\int_{-b}^{-1}\left|\left(\Phi(x+z)-\Phi(x)\right)\frac{{{C^{j}_{-}}w(-{\beta^{j}_{-}}z)}}{{{{(-z)}^{1+{\alpha^{j}_{-}}}}}}\right|\mathrm{d}z$
$\displaystyle\leq
2\Phi(0^{+})\int_{-b}^{-1}\frac{C^{j}_{-}w(-{\beta^{j}_{-}}z)}{{{{(-z)}^{1+{\alpha^{j}_{-}}}}}}\mathrm{d}z=\Phi(0^{+})\int_{b>|z|\geq
1}\nu_{\rm sym}(j,\mathrm{d}z),$

in which $\nu_{\rm sym}(j,\mathrm{d}z)$ is the “symmetrized” version of the
original Lévy measure. Moreover, we have

$\displaystyle\left|I^{(4)}_{v}\right|$
$\displaystyle\leq\int_{1}^{b}\left|\left(\Phi(x+z)-\Phi(x)\right)\frac{{{C^{j}_{+}}w({\beta^{j}_{+}}z)}}{{{{z}^{1+{\alpha^{j}_{+}}}}}}\right|\mathrm{d}z\leq\Phi(0^{+})\int_{1\leq|z|<b}\nu_{\rm
sym}(j,\mathrm{d}z).$

According to Remark 2.2, the right hand sides of the above two inequalities
are absolutely convergent as $b\rightarrow\infty$. Note also that the
symmetric Lévy measure $\nu_{\rm sym}(j,\mathrm{d}z)$ is given by substituting
$C^{j}_{+}=C^{j}_{-}$, $\alpha^{j}_{+}=\alpha^{j}_{-}$ and
$\beta^{j}_{+}=\beta^{j}_{-}$ in (2.8). The proof will be complete if we show
that the integrals $I^{(2)}_{v}$ and $I^{(3)}_{v}$ are also absolutely
convergent. In this respect, we will simultaneously prove them for the case
$\alpha_{\pm}<2$, which consists of the following three different subcases:

  * (1)

Finite activity tempered stable Lévy measures ($\alpha_{\pm}<0$),

  * (2)

Infinite activity but finite variation tempered stable Lévy measures
($0\leq\alpha_{\pm}<1$),

  * (3)

Infinite activity and infinite variation tempered stable Lévy measures
($1\leq\alpha_{\pm}<2$).

Based on Definition 2.1 and also by applying Theorem 4.6, Theorem 4.7, Lemma
4.8 and also Taylor’s theorem, the following inequalities are obtained for the
integrals $I^{(2)}_{v}$ and $I^{(3)}_{v}$:

$\displaystyle\left|I^{(2)}_{v}\right|$
$\displaystyle\leq\int^{-\varepsilon}_{-1}\left|\left(\Phi(x+z)-\Phi(x)-z\frac{\partial}{\partial
x}\Phi\right)\frac{{{C^{j}_{-}}w({-\beta^{j}_{-}}z)}}{{{{(-z)}^{1+{\alpha^{j}_{-}}}}}}\right|\mathrm{d}z$
$\displaystyle=\int^{-\varepsilon}_{-1}\left|\left(\frac{z^{2}}{2!}\frac{\partial^{2}}{\partial
x^{2}}\Phi(x)+R_{2}(z)\right)\frac{{{C^{j}_{-}}w({-\beta^{j}_{-}}z)}}{{{{(-z)}^{1+{\alpha^{j}_{-}}}}}}\right|\mathrm{d}z$
$\displaystyle\leq\frac{M_{2}}{2!}\int^{-\varepsilon}_{-1}z^{2}\frac{{{C^{j}_{-}}w({-\beta^{j}_{-}}z)}}{{{{(-z)}^{1+{\alpha^{j}_{-}}}}}}\mathrm{d}z+\frac{M_{3}}{3!}\int^{-\varepsilon}_{-1}|z|^{3}\frac{{{C^{j}_{-}}w({-\beta^{j}_{-}}z)}}{{{{(-z)}^{1+{\alpha^{j}_{-}}}}}}\mathrm{d}z$
$\displaystyle\leq\frac{M_{2}}{4}\int_{\varepsilon\leq|z|<1}z^{2}\nu_{\rm
sym}(j,\mathrm{d}z)+\frac{M_{3}}{12}\int_{\varepsilon\leq|z|<1}|z|^{3}\nu_{\rm
sym}(j,\mathrm{d}z)$
$\displaystyle\leq\frac{1}{4}\max\\{M_{2},M_{3}\\}\int_{\varepsilon\leq|z|<1}z^{2}\nu_{\rm
sym}(j,\mathrm{d}z),$

and

$\displaystyle\left|I^{(3)}_{v}\right|$
$\displaystyle\leq\int^{1}_{\varepsilon}\left|\left(\Phi(x+z)-\Phi(x)-z\frac{\partial}{\partial
x}\Phi\right)\frac{{{C^{j}_{+}}w({\beta^{j}_{+}}z)}}{{{{z}^{1+{\alpha^{j}_{+}}}}}}\right|\mathrm{d}z$
$\displaystyle=\int^{1}_{\varepsilon}\left|\left(\frac{z^{2}}{2!}\frac{\partial^{2}}{\partial
x^{2}}\Phi(x)+R_{2}(z)\right)\frac{{{C^{j}_{+}}w({\beta^{j}_{+}}z)}}{{{{z}^{1+{\alpha^{j}_{+}}}}}}\right|\mathrm{d}z$
$\displaystyle\leq\frac{M_{2}}{2!}\int_{\varepsilon}^{1}z^{2}\frac{{{C^{j}_{+}}w({\beta^{j}_{+}}z)}}{{{{z}^{1+{\alpha^{j}_{+}}}}}}\mathrm{d}z+\frac{M_{3}}{3!}\int_{\varepsilon}^{1}|z|^{3}\frac{{{C^{j}_{+}}w({\beta^{j}_{+}}z)}}{{{{z}^{1+{\alpha^{j}_{+}}}}}}\mathrm{d}z$
$\displaystyle\leq\frac{M_{2}}{4}\int_{\varepsilon\leq|z|<1}z^{2}\nu_{\rm
sym}(j,\mathrm{d}z)+\frac{M_{3}}{12}\int_{\varepsilon\leq|z|<1}|z|^{3}\nu_{\rm
sym}(j,\mathrm{d}z)$
$\displaystyle\leq\frac{1}{4}\max\\{M_{2},M_{3}\\}\int_{\varepsilon\leq|z|<1}z^{2}\nu_{\rm
sym}(j,\mathrm{d}z),$

in which $R_{k}(z)$ is the mean-value form of the remainder term in the
Taylor’s theorem and given by:

$R_{k}(z)=\frac{f^{(k+1)}(\xi_{x})}{(k+1)!}z^{k+1},$ (4.17)

for some real number $\xi_{x}$ between $x$ and $x+z$.

By the comparison test for improper integrals and according to the definition
of a Lévy measure (see Subsection 2.1), the integrals $I_{v}^{(2)}$ and
$I_{v}^{(3)}$ are absolutely convergent due to the existence of upper integral
bounds as $\varepsilon\rightarrow 0^{+}$. Note that
$|R_{k}(z)|\leq\frac{||f^{(k+1)}(\cdot)||}{(k+1)!}|z|^{k+1}$ and so the proof
is complete. ∎

Theorem 4.9 shows that collocation based on infinitely differentiable radial
basis functions such as Gaussian bases is a reliable method to approximate the
solution of PIDEs arising from models including tempered stable Lévy
processes. I the next section, we assess this in practice for different
standard test problems.

##  5 Numerical Experiments

In this section, we use numerical experiments to showcase the use of the
proposed meshfree methodology to find the default probability of firms whose
asset price model is a Markov-modulated Lévy process with synchronous jumps.
Uniform distribution of points in interval $[x_{\rm min},x_{\rm max}]=[-8,8]$
is used for the collocation procedure [13]. Although, we analyze the
methodology with different families of RBFs, we choose to report the results
based on Gaussian RBFs due to the fact that the solution of a PIDEs based on
the jump-diffusion model inherently behaves like a Gaussian distribution (see
[56] for more details). The first test problem is pertained to three French
firms calibrated and fitted by Markov modulated Lévy processes with
synchronous jumps in [11] for a period of 10 years from the first of January
2004 to the twenty forth of January 2014, whereas the two last test problems
are designed by authors for CGMY Lévy process and Kobol lévy process that are
analyzed with the proposed method in this paper.

Test Problem I (VG Model) In the first experiment, our computations are based
on parameters obtained in [11] for three French companies: Axa,
STMicroelectonics and société Générale with the fitted parameters reported in
Table 2. Using econometric adjustment in [11], a regime switching process with
two regimes ($H=2$) based on the VG process with synchronous jumps is used to
model the dynamics of the firms’s asset value process.

Table 2: Parameters of the VG model for Axa, STMicroelectonics and Société Générale firms. VG Parameters | Axa |  | STM |  | Soc Gen
---|---|---|---|---|---
| State 1 | State 2 |  | State 1 | State 2 |  | State 1 | State 2
$\sigma_{j}$ | 0.4460 | 0.1234 |  | 0.2495 | 0.1534 |  | 0.3227 | 0.1675
$\theta_{j}$ | -0.1421 | 0.0196 |  | -0.1135 | 0.0043 |  | -0.1576 | 0.0254
$\kappa_{j}$ | 0.0236 | 0.0011 |  | 0.0374 | 0.0015 |  | 0.0306 | 0.0028
$p_{j}$ | 0.8755 | 0.9664 |  | 0.6822 | 0.9386 |  | 0.8083 | 0.9549
$\eta_{i,j}$ | 0.0160 | -0.0092 |  | 0.0078 | -0.0074 |  | 0.0132 | -0.0117

For this problem, the probability of default surfaces for the
STMicrolectronics firm are illustrated in two different regimes in Figure 4.
Besides, we compare the results of meshfree collocation method of this paper
with the result of FFT method proposed by Hainaut et al. in [11] in Table 3.
It can be seen that we obtain better results as the number of spatial nodes
increases. Nevertheless the computational time will naturally increase as the
size of the system becomes larger. The probability of bankruptcy for the three
companies is depicted in Figure 2 for the time horizon from 1 to 10 years,
where firms can default at the end of each financial year.

Figure 1:  Figure 2: Default Probability, VG models Table 3: The error incurred in estimating the probability of defaulat for Société Générale. $N_{s}$ | Relative Error in State 1 | Relative Error in State 2 | CPU-Time (s)
---|---|---|---
64 | 0.638564 | 0.67354 | 0.814
128 | 0.012792 | 0.01321 | 0.955
256 | 0.053485 | 0.05090 | 1.812
512 | 0.001012 | 0.00103 | 2.336

(a)

(b)

Figure 3: Accuracy properties of numerical method for Axa with VG models, a)
Impact of space steps b) Impact of time steps

(a)

(b)

Figure 4: Probability of default for STEM: (a) the first regime (b) the second
regime. Figure 5: Default Probability for the CGMY model. Figure 6:
Probability of default based on the CGMY Process for state 1.

In order to study the numerical error of the meshfree RBF collocation method
for evaluating the probability of default with respect to discretization
parameters of $\Delta x$ and $\Delta\tau$, we calculate the log relative error
obtained by

$\|E\|_{\text{LRE}}=\log\Big{|}\frac{\text{PD}_{\text{RBF}}-\text{PD}_{\text{FFT}}}{\text{PD}_{\text{FFT}}}\Big{|}$

In Figure 3, we use this error metric to assess the rate of convergence of our
RBF collocation method for evaluating the probability of default for the Axa
case. It is evident from this figure that the rate of convergence of the
numerical scheme is almost linear in time and super-linear in space. In
general, it seems that the integral term related to the synchronous jump term
(3.4) the accuracy, and deteriorates the rate of convergence in time and
specially space.

Based on the fact that there is no real world example of calibrating the asset
of a firm in tempered stable Lévy model such as CGMY and Kobol models in the
form of a regime-switching synchronous jump Lévy process, so we have designed
two artificial test cases based on these two processes similar to the example
designed in [13].

Test Problem II (CGMY Model) For this example, we assume that the value of a
hypothetical firm follows a CGMY Lévy process in the context of regime-
switching synchronise-jump dynamics containing 5 different regimes. Let us
also assume that the following generator matrix, $Q$, for the transition
probability matrix is given by

$Q=\begin{bmatrix}-1&0.25&0.25&0.25&0.25\\\ 0.25&-1&0.25&0.25&0.25\\\
0.25&0.25&-1&0.25&0.25\\\ 0.25&0.25&0.25&-1&0.25\\\ 0.25&0.25&0.25&0.25&-1\\\
\end{bmatrix}$ (5.1)

Data related to this model for five different states is provided in Table 4.

Table 4: Parameters of the CGMY model. CGMY Parameters | State 1 | State 2 | State 3 | State 4 | State 5
---|---|---|---|---|---
$C$ | 0.1 | 0.3 | 0.5 | 0.7 | 0.9
$G$ | 2.0 | 4.0 | 6.0 | 8.0 | 10.0
$M$ | 1.0 | 3.0 | 5.0 | 7.0 | 9.0
$Y$ | 0.11 | 0.22 | 0.33 | 0.44 | 0.55
$\eta_{i,j}$ | 0.01 | -0.02 | 0.03 | -0.04 | 0.05

We have used a computational domain of the form, $[x_{\rm min},x_{\rm
max}]=[-10,10]$ for our collocation method and also for truncating the
integration domain. The approximated probability of default for this case and
for time periods from 1 to 10 years is depicted in Figure (5). It must be
noted that besides the Gaussian basis function, the cubic radial basis
function, $\phi(r)=r^{3}$, works perfectly well and produces reliable and
stable solutions. Figure (6) portraits the portability of default surface in
this problem for the first regime.

Test Problem III (KoBoL Model) This test case is dedicated to finding the
probability of default, when the dynamics of the firm’s asset follows a
regime-switching Lévy process with three different regimes based on KoBoL jump
process and synchronous jumps effect. The generator matrix, $Q$ of the
transition probability matrix is given by

$Q=\begin{bmatrix}-1&0.3&0.7\\\ 0.3&-1&0.7\\\ 0.7&0.3&-1\end{bmatrix},$ (5.2)

with model parameters provided in Table (5).

Table 5: Parameters of the KoBol model. KoBol Parameters | State 1 | State 2 | State 3
---|---|---|---
$C$ | 0.1 | 0.11 | 0.13
$Y$ | 0.9 | 1.2 | 1.8
$P$ | 0.2 | 0.4 | 0.8
$q$ | 0.3 | 0.5 | 0.7
$\lambda$ | 3.0 | 2 | 2.5
$\eta_{i,j}$ | 0.04 | -0.01 | 0.02

Computational domain is $[x_{\rm min},x_{\rm max}]=[-10,10]$ for collocation
points. The probability of default of the firm for three different regimes is
illustrated in Figure (7) in the time horizon from 1 to 10 years when the
default is happening just at the maturity.

Figure 7: Default probability for KoBol model.

Moreover, Figure (8(c)) depicts the probability of default surface in this
case with three different regimes with the KoBol jump model over the space-
time domain of $S\times t\in[e^{-10},e^{10}]\times[0,10]$.

(a)

(b)

(c)

Figure 8: Probability of default in KoBol Model: (a) First regime; (b) Second
regime; (c) Third regime

.

The stable behavior of these surfaces confirm that the proposed method is
reliable and useful for problems with tempered stable jump terms as well as
problems with more economic regimes.

##  6 Conclusion

In this investigation, a radial basis function collocation method is proposed
to evaluate the default probability of a firm where the value of the firm is
satisfying the newly proposed synchronous-jump regime-switching model. We
explained the detailed methodology of the proposed method including deriving a
semi-discrete system from discretizing the time derivative with finite-
differences and meshfree collocation on the spatial direction. Generalizing
the proposed synchronous-jump regime-switching model to benefit from tempered
stable processes, we theoretically proved that the proposed numerical method
is stable and inherently de-singularizes the problem. The numerical
experiments confirm the efficiency and accuracy of the method. This method can
be utilized for other structural models based on the synchronous-jump regime-
switching model both for pricing financial derivatives and also credit risk
computations.

## References

  * [1] A. J. McNeil, R. Frey, P. Embrechts, Quantitative Risk Management: Concepts, Techniques and Tools, Princeton University Press, 2015.
  * [2] H. G. Fong, The Credit Market Handbook: Advanced Modeling Issues, Vol. 340, John Wiley & Sons, 2006.
  * [3] D. Duffie, K. J. Singleton, Credit risk: pricing, measurement, and management, Princeton university press, 2012.
  * [4] F. Modigliani, M. H. Miller, The cost of capital, corporation finance and the theory of investment, The American economic review 48 (3) (1958) 261–297.
  * [5] F. Black, M. Scholes, The pricing of options and corporate liabilities, Journal of Political Economy 81 (3) (1973) 637–654.
  * [6] R. C. Merton, On the pricing of corporate debt: The risk structure of interest rates, The Journal of finance 29 (2) (1974) 449–470.
  * [7] S. Sundaresan, A review of merton’s model of the firm’s capital structure with its wide applications, Annu. Rev. Financ. Econ. 5 (1) (2013) 21–41.
  * [8] F. Fiorani, E. Luciano, P. Semeraro, Single and joint default in a structural model with purely discontinuous asset prices, Quantitative Finance 10 (3) (2010) 249–263.
  * [9] T. K. Siu, C. Erlwein, R. S. Mamon, The pricing of credit default swaps under a markov-modulated merton’s structural model, North American Actuarial Journal 12 (1) (2008) 18–46.
  * [10] T. S. R.J. Elliott, L. Chan, Option pricing and Esscher transform under regime switching, Annals of Finance (2005) 423–432.
  * [11] D. Hainaut, D. B. Colwell, A structural model for credit risk with switching processes and synchronous jumps, The European Journal of Finance 22 (11) (2016) 1040–1062.
  * [12] K. Chourdakis, Switching lévy models in continuous time: Finite distributions and option pricing, University of Essex, Centre for Computational Finance and Economic Agents (CCFEA) Working Paper.
  * [13] A. F. Bastani, Z. Ahmadi, D. Damircheli, A radial basis collocation method for pricing american options under regime-switching jump-diffusion models, Applied Numerical Mathematics 65 (2013) 79–90.
  * [14] Y. Lee, Financial options pricing with regime-switching jump-diffusions, Computers & Mathematics with Applications 68 (3) (2014) 392–404.
  * [15] N. Rambeerich, A. A. Pantelous, A high order finite element scheme for pricing options under regime switching jump diffusion processes, Journal of Computational and Applied Mathematics 300 (2016) 83–96.
  * [16] P. A. A. Tonfack, Modeling in finance and insurance with Levy-Itö driven dynamic processes under semi markov-type switching regimes and time domains, Ph.D. thesis, Ph.D. Thesis, University of South Florida (2017).
  * [17] D. Hainaut, D. B. Colwell, A structural model for credit risk with switching processes and synchronous jumps, The European Journal of Finance 22 (11) (2016) 1040–1062.
  * [18] K. R. Jackson, S. Jaimungal, V. Surkov, Option pricing with regime switching lévy processes using fourier space time stepping, in: Proc. 4th IASTED Intern. Conf. Financial Engin. Applic, 2007, pp. 92–97.
  * [19] K. R. Jackson, S. Jaimungal, V. Surkov, Fourier space time-stepping for option pricing with lévy models, Journal of Computational Finance 12 (2) (2008) 1–29.
  * [20] K.-I. Sato, Lévy Processes and Infinitely Divisible Distributions, Cambridge University Press, 1999.
  * [21] U. Küchler, S. Tappe, Tempered stable distributions and processes, Stochastic Processes and their Applications 123 (12) (2013) 4256–4293.
  * [22] P. Carr, H. Geman, D. B. Madan, M. Yor, The fine structure of asset returns: An empirical investigation, The Journal of Business 75 (2) (2002) 305–332.
  * [23] D. B. Madan, E. Seneta, The variance gamma (vg) model for share market returns, Journal of business (1990) 511–524.
  * [24] D. B. Madan, Purely discontinuous asset price processes, Option pricing, Interest rates and risk management (2001) 105–153.
  * [25] A. Lee, Y. Lee, Comparison of numerical methods for option pricing under the cgmy model, Journal of the Chungcheong Mathematical Society 29 (3) (2016) 503–508.
  * [26] D. Revuz, M. Yor, Continuous Martingales and Brownian Motion, Vol. 293, Springer Science & Business Media, 2013.
  * [27] M. L. Bianchi, S. T. Rachev, Y. S. Kim, F. J. Fabozzi, Tempered stable distributions and processes in finance: numerical analysis, in: Mathematical and statistical methods for actuarial sciences and finance, Springer, 2010, pp. 33–42.
  * [28] S. T. Rachev, Y. S. Kim, M. L. Bianchi, F. J. Fabozzi, Financial models with Lévy processes and volatility clustering, Vol. 187, John Wiley & Sons, 2011\.
  * [29] P. Tankov, Financial modelling with jump processes, Chapman and Hall/CRC, 2003.
  * [30] U. Küchler, S. Tappe, Bilateral gamma distributions and processes in financial mathematics, Stochastic Processes and their Applications 118 (2) (2008) 261–283.
  * [31] A. Almendral, C. W. Oosterlee, Accurate evaluation of european and american options under the cgmy process, SIAM Journal on Scientific Computing 29 (1) (2007) 93–117.
  * [32] A. Hirsa, Computational Methods in Finance, CRC Press, 2012.
  * [33] R. Cont, E. Voltchkova, A finite difference scheme for option pricing in jump diffusion and exponential lévy models, SIAM Journal on Numerical Analysis 43 (4) (2005) 1596–1626.
  * [34] R. Brummelhuis, R. T. Chan, A radial basis function scheme for option pricing in exponential levy models, Applied Mathematical Finance 21 (3) (2014) 238–269.
  * [35] A. Itkin, Pricing Derivatives under Lévy Models, Vol. 12, Springer, 2017.
  * [36] D. Damircheli, M. Bhatia, Solution approaches and sensitivity analysis of variational inequalities, in: AIAA Scitech 2019 Forum, 2019, p. 0977.
  * [37] C. La Chioma, Integro-differential problems arising in pricing derivatives in jump-diffusion markets, Ph.D. thesis, Ph. D. Thesis.(Roma, 2003-2004) (2003).
  * [38] V. Surkov, Option pricing using fourier space time-stepping framework, Ph.D. thesis, Ph.D. Thesis, University of Toronto (2009).
  * [39] G. E. Fasshauer, Meshfree approximation methods with MATLAB, Vol. 6, World Scientific, 2007.
  * [40] K. Kormann, E. Larsson, A Galerkin radial basis function method for the Schrödinger equation, SIAM Journal on Scientific Computing 35 (6) (2013) 2832–2855.
  * [41] S. A. Sarra, A numerical study of the accuracy and stability of symmetric and asymmetric rbf collocation methods for hyperbolic PDEs, Numerical Methods for Partial Differential Equations: An International Journal 24 (2) (2008) 670–686.
  * [42] P. Brémaud, Markov Chains: Gibbs Fields, Monte Carlo Simulation and Queues, Springer-Verlag, 1999.
  * [43] A. Papantaleon, An introduction to lévy processes with applications in finance, University of Freiburg.
  * [44] R. C. Merton, Option pricing when underlying stock returns are discontinuous, Journal of financial economics 3 (1-2) (1976) 125–144.
  * [45] S. G. Kou, A jump-diffusion model for option pricing, Management science 48 (8) (2002) 1086–1101.
  * [46] S. I. Boyarchenko, S. Z. Levendorskiǐ, Option pricing for truncated lévy processes, International Journal of Theoretical and Applied Finance 3 (03) (2000) 549–552.
  * [47] I. Koponen, Analytic approach to the problem of convergence of truncated lévy flights towards the gaussian stochastic process, Physical Review E 52 (1) (1995) 1197.
  * [48] X. Guo, R. A. Jarrow, Y. Zeng, Modeling the recovery rate in a reduced form model, Mathematical Finance: An International Journal of Mathematics, Statistics and Financial Economics 19 (1) (2009) 73–97.
  * [49] F. A. Bornemann, Adaptive multilevel discretization in time and space for parabolic partial differential equations, Bericht TR 89-07, Konrad-Zuse-Zentrum Berlin.
  * [50] G. E. Fasshauer, Newton iteration with multiquadrics for the solution of nonlinear pdes, Computers & Mathematics with Applications 43 (3-5) (2002) 423–438.
  * [51] A. F. Bastani, M. V. Dastgerdi, A. Mighani, On multilevel rbf collocation to solve nonlinear pdes arising from endogenous stochastic volatility models, Communications in Nonlinear Science and Numerical Simulation 59 (2018) 88–104.
  * [52] S. Bochner, Harmonic Analysis and the Theory of Probability, University of California Press, 1955.
  * [53] J. Buescu, A. C. Paixão, On differentiability and analyticity of positive definite functions, Journal of mathematical analysis and applications 375 (1) (2011) 336–341.
  * [54] E. Massa, A. Peron, A. Piantella, Estimates on the derivatives and analyticity of positive definite functions on $\mathbb{R}^{m}$, Analysis Mathematica 43 (1) (2017) 89–98.
  * [55] R. Wrede, M. Spiegel, Schaum’s Outline of Advanced Calculus, McGraw-Hill Education, 2002.
  * [56] S.-M.-M. Kazemi, M. Dehghan, A. Foroush Bastani, On a new family of radial basis functions: Mathematical analysis and applications to option pricing, Journal of Computational and Applied Mathematics 328 (2018) 75–100.