{'fraction_non_alphanumeric': 0.037868162692847124, 'fraction_numerical': 0.002428762015530394, 'mean_word_length': 4.4308006687720605, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Meshfree Collocation to Compute the Probability of Default under a Regime-Switching Synchronous-Jump Tempered Stable L\'{e}vy Model

In the paper [Hainaut, D. and Colwell, D.B., {\rm 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\'{e}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.


. 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 accords , 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 [ ]. They also provide useful tools for market participants to identify, measure, monitor and control their exposure to credit risk [ ].
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 [ ]. 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 [ ] and option pricing theory of Black and Scholes [ ] and Merton [ ] and rely heavily on diffusion processes to model the evolution of the firm's asset value process [ ].
There have been some efforts in the literature to extend the Merton's structural model to more complex dynamical processes such as jump-diffusion [ ] or regime-switching diffusion processes [ ].
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. [ , ] 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 [ ] for a complete account) with some option pricing applications (see e.g. [ , , ]). 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. [ ]).
In a recent contribution, Hainaut and Colwell [ ] 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. [ , ]) 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 [ ] 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 [ ], where the associated Lévy processes are called "truncated Lévy flights" and have been generalized by several authors [ ]. 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 [ , ], bilateral Gamma distributions [ ] and CGMY distributions [ ]. Properties of tempered stable distributions have been investigated, e.g., in [ , , , ]. For financial modelling, they have been applied, e.g. in [ , , , ], see also the recent book [ ].
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. [ , , , , , , , ]), 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 [ , , ] have approximated the above-mentioned PIDEs with a new singularly perturbed PIDE( ) 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 [ ]. 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 [ ]. 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 [ ]. 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. [ , ] 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 [ ] 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 [ ] 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 , 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 , 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 present the numerical method based on the proposed radial basis function collocation. In Section , 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.

. 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, α t , with values from the set H = {1, 2, · · · , H}, each representing an economic state or regime. Let Q = (q ij ) H×H denote the generator matrix of α t where the off-diagonal entries are non-negative and H j=1 q ij = 0, for all i ∈ H (see e.g. [ ]). Discretizing the time variable with a sufficiently small step-size, ∆t, we could show that q i,j ∆t is the probability to switch from state i to state j, for i = j. We also define the transition probability matrix of the process, denoted by P (t, s) in the form 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.

. . 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 in which X t is a Markov-modulated Lévy process with X 0 = 0. More precisely, X t is a piecewise Lévy process depending upon the state, α t of an economy with simultaneous jumps at regime switching times. Dynamics of the process X t could be written as where X t has jumps represented by J α − t ,αt when regime switches at time t. The processes X j t , j = 1, · · · H are independent Lévy processes with a Lévy-Itô decomposition of the form 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 ν(j, ·) (different from J α − t ,αt ), the Lévy measure of X t in state j. The triplet (µ j , σ j , ν(j, z)) uniquely determines the characteristic function of X j t by the Lévy-Khintchine theorem (see e.g. [ , ]) stating that in which the characteristic exponent of the Lévy process is given by In the above expressions, µ j is the drift, σ j > 0 is the diffusion and ν(j, ·) is the Lévy measure on R satisfying ν(j, {0}) = 0 and R (1 ∧ |z| 2 )ν(j, dz) < ∞ (see e.g. [ ] for more details on Lévy processes).
The finite-activity Lévy density of synchronous jumps J i,j is given by (see [ ]) which can take positive, negative or real values with a characteristic function of the form In the next subsection, we introduce the specific Lévy process used in this paper along with its Lévy measure.

. . Generalized Tempered Stable Lévy Processes
Consider a Lévy process, X j t , specified by it's characteristic triplet, (µ j , σ j , ν(j, z)). The process is called of finite activity when ν(j, R) = R ν(j, dz) < ∞ and of infinite activity otherwise (see e.g. [ ]). We know also that a Lévy process can be decomposed in the form where J j X describes a jump process that its jumps can be finite or infinite. Let the family of functions, W, to contain those w : (0, ∞) → [0, ∞) which are continuous and also satisfy the following conditions: • They are decreasing, Definition . . A Lévy process is of generalized tempered-stable (GTS) type if its Lévy measure, ν(j, ·), could be expressed as in which w ∈ W and the parameters α j ± < 2 and C j ± , β j ± > 0 are chosen such that |z| n C j ± w(β j ± z) is bounded on R for each n ∈ N.

Indeed, the relation ( . ) is a necessary and sufficient condition for moments of a GTS Lévy process to exist [ ].
Remark . . The following particular cases are known in the literature: [ , ]; .

Different types of Lévy jump-diffusion processes and their Lévy-Khintchin representation
If a Lévy jump-diffusion process is of finite activity, then it has also jumps of finite variation. .

. 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 (Ω, F, Q) equipped with some filtration, 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 Q is taken to be a risk neutral measure but the following results are also valid under the real historical measure, 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: and the probability of default is denoted by In case of bankruptcy, a fraction R(α T ), called the recovery rate, of asset value is repaid to debtholders.
In this work, as in Guo et al. [ ], the debt recovers at a different magnitude depending here upon the economic regime, α(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 ( . ) We will present the approach used to compute the probabilities of default and the defaultable bond prices in the reminder.

. 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: . Using Itô's formula for general semi-martingales (or infinitesimal generator) [ ]; . Using the Lévy-Khintchine formula and the tempered Fourier transform (see Proposition The second approach is more useful in characterizing some intrinsic properties of the integrodifferential operator of the problem which will be needed in later parts of the paper. The probability of default of a firm in state α t and time t under the measure P is defined as in which D (t, T, x, j) 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 [ ]) that D (t, T, x, j) for j = 1, 2, · · · , H satisfy the following system of partial integro-differential equations in which L is the infinitesimal generator of default probabilities of the form (see e.g. [ , ]) Due to the relation H i=1 q j,i = 0, the generator can be written as In order to make the presentation and implementation of the numerical meshfree scheme in the next section simple, we transform the PIDEs ( . ) forward in time. To this end and using a common change of variables of the form t = T − τ and thus u(τ, x, j) = D (t, T, x, j) where T is removed from u for simplicity of notation, the problem ( . ) will be transformed into Developing a stable numerical scheme to solve Eq. ( . ) 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.

. De-Singularized Meshfree Approximation
The important point to notice in confronting with the infinite activity processes is the singularity of the Lévy measure ν(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 α j ± > 0) approaches infinity much faster than the variance gamma (VG) process [ ]. For this reason we should have the highest possible order of approximation for the integrand near zero.
Indeed, as we move in Table 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. [ , , , , , , ]).
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 , PIDE( ) is used to refer to this derived PIDE. When goes to 0, it is proved that the solution of PIDE( ) will converge to the solution of original PIDE ([ , , ]). Applying this approach to the system ( . ) will lead us to the following PIDE( ) where the factors ω j (ε) and σ j (ε) have been introduced in [ ] (Chapter ).
• 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 [ ]).
• Last but not least, we will mention the method used by

. . 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 ( . ) which could be represented as in which 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. [ ]). We follow this tradition by first discretizing the equation ( . ) 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 [ ] and [ ] for more details).

. Meshfree RBF Collocation
In order to give a rough idea of the radial basis function collocation, we first need the following concepts: Note that the norms will be the usual Euclidean -norm, defined by x 2 := √ x T 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 (τ, x) onto a finite dimensional space as a linear combination of radial basis function of the form in which, x j for j = 1, · · · , N s is a set of scattered data points in the sub domain I and Υ So we have It is worth mentioning that mega matrix Φ L has a block structural where off-diagonal blocks entangles data from different regimes (see for more details [ ]). Let us recall some fundamental properties of RBFs before presenting the core theorem of this section.

. . 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 . . (Bochner [ ]) A radially symmetric multivariate function
The function Φ is called strictly positive definite on R d , if the quadratic form in the left hand side of ( . ) is zero only for c = 0. We also call a univariate function φ : [0, ∞) → R strictly positive definite on R d (abbreviated as φ ∈ SPD d ), if the corresponding radially symmetric multivariate function defined by Φ(x) := φ( x ), x ∈ R d is strictly positive definite.
In order to prove our main result, we need the following properties of SPD d functions which play a key role in our development.

Theorem . . (Buescu and Paixao [ ] and Massa and et al. [ ])
Let Φ : R m → R be a positive definite function and suppose Φ is of class C 2n (resp. C ∞ ) in some neighborhood of the origin for some positive integer n. Then • In the special case where m = 1, we have Φ ∈ C 2n (R) and for all integers m 1 , m 2 with 0 ≤ m i ≤ n, i = 1, 2 and every x ∈ R we have

Remark . . Based on the above results, we have the following:
(i) Every SPD d function, Φ, is bounded above (see the Theorem . ).

(iii) Based on the Theorem . , if we define
Theorem . . Let Φ : R → R be a real-valued strictly positive definite function and assume that Φ is of class C 2n in some neighborhood of the origin for some positive integer n ≥ 2. Then, we have Φ ∈ C 2n (R). Also for any tempered stable Lévy measure ν(i, dz) with α j ± < 2 defined in . , the integral is finite for all x ∈ R and without singularity at the origin.
Proof. For arbitrary constants 0 < ε < 1 and b 1, the improper integral I v (·) could be written as According to the comparison test for improper integrals (see Chapter in [ ]), it is easy to show that the improper integrals of the first kind (i.e. I (1) v and I (4) v ) are absolutely convergent as b tends to infinity for arbitrary α j ± < 2, because we can obtain the following upper bounds for the integrals I (1) v and I (4) v using Theorem . and Definition . : in which ν sym (j, dz) is the "symmetrized" version of the original Lévy measure. Moreover, we have According to Remark . , the right hand sides of the above two inequalities are absolutely convergent as b → ∞. Note also that the symmetric Lévy measure ν sym (j, dz) is given by substituting . The proof will be complete if we show that the integrals I v are also absolutely convergent. In this respect, we will simultaneously prove them for the case α ± < 2, which consists of the following three different subcases: Theorem . 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.

. 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 Markovmodulated Lévy process with synchronous jumps. Uniform distribution of points in interval [x min , x max ] = [−8, 8] is used for the collocation procedure [ ]. 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 [ ] 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 [ ] for a period of years from the first of January to the twenty forth of January , 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 [ ] for three French companies: Axa, STMicroelectonics and société Générale with the fitted parameters reported in Table . Using econometric adjustment in [ ], 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.
For this problem, the probability of default surfaces for the STMicrolectronics firm are illustrated in two different regimes in Figure . Besides, we compare the results of meshfree collocation method of this paper with the result of FFT method proposed by Hainaut et al. in [ ] in Table . 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 for the time horizon from to years, where firms can default at the end of each financial year. 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 ∆x and ∆τ , we calculate In Figure , 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 ( . ) 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 [ ]. 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 different regimes. Let us also assume that the following generator matrix, Q, for the transition probability matrix is given by Data related to this model for five different states is provided in Table . We have used a com-  Figure ( ). It must be noted that besides the Gaussian basis function, the cubic radial basis function, φ(r) = r 3 , works perfectly well and produces reliable and stable solutions. Figure   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 with model parameters provided in Table ( ). Computational domain is [x min , x max ] = [−10, 10] for collocation points. The probability of default of the firm for three different regimes is illustrated in Figure ( ) in the time horizon from to years when the default is happening just at the maturity. Moreover, Figure ( ) depicts the probability of default surface in this case with three different  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.

. 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 synchronousjump 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.