PDE models for the valuation of a non callable defaultable coupon bond under an extended JDCEV model

TL;DR

This paper develops PDE-based numerical methods for valuing non-callable, defaultable bonds under an extended JDCEV model that unifies credit and equity risks, with validation against Monte Carlo and asymptotic methods.

Contribution

It introduces a PDE framework with specific numerical schemes for bond valuation under the JDCEV model, including existence proofs and real-market numerical validation.

Findings

Numerical schemes accurately price bonds consistent with market data.

The PDE approach aligns well with Monte Carlo and asymptotic methods.

The model captures complex credit and equity risk interactions.

Abstract

We consider a two-factor model for the valuation of a non callable defaultable bond which pays coupons at certain given dates. The model under consideration is the Jump to Default Constant Elasticity of Variance (JDCEV) model. The JDCEV model is an improvement of the reduced form approach, which unifies credit and equity models into a single framework allowing for stochastic and possible negative interest rates. From the mathematical point of view, the valuation involves two partial differential equation (PDE) problems for each coupon. First, we obtain the existence of solution for these PDE problems. In order to solve them, we propose appropriate numerical schemes based on a Crank-Nicolson semi-Lagrangian method for time discretization combined with bi-quadratic Lagrange finite elements for space discretization. Once the numerical solutions of the PDEs are obtained, a post-processing procedure is carried out in order to achieve the value of the bond. This post-processing includes the computation of an integral term which is approximated by using the composite trapezoidal rule. Finally, we present some numerical results for real market bonds issued by different firms in order to illustrate the proper behaviour of the numerical schemes. Moreover, we obtain an agreement between the numerical results from the PDE approach and those ones obtained by applying a Monte Carlo technique and an asymptotic approximation method.

Figures

Tables

Table 1. Table 1: Different finite element meshes (number of elements and nodes).

	Number of elements	Number of nodes
Mesh 4	16	81
Mesh 8	64	289
Mesh 16	256	1089
Mesh 32	1024	4225

Table 2. Table 2: Parameters of the model for the UBS bond.

Parameters of the defaultable stock price model

$a_1=0.0337851$

$a_2=0.0523625$

$b_1=0.0026639$

$b_2=0.0027968$

$c=0.0435673$

$\beta =-0.268496$

Parameters of the interest rate model

$\kappa =0.04520533766268042$

$\delta =0.02146900332086033$

$\theta =0.10334921942765922$

Correlation coefficient

$\rho =0.0$

Initial conditions

$S_0=1.0$

$r_0=-0.009159871729892612$

Table 3. Table 3: Market and model values of the ZCB.

Maturity (years)	Market ZCB	Model ZCB
1	1.00229	1.006751
2	1.00372	1.009062
3	1.00333	1.007495
4	1.00099	1.002601
5	0.995825	0.994902
6	0.987805	0.984889
7	0.976833	0.973024
8	0.963223	0.959738
9	0.947687	0.945429
10	0.932845	0.930463

Table 4. Table 4: Value of the UBS bond for different meshes and time steps.

Time steps per year	Mesh 4	Mesh 8	Mesh 16	Mesh 32
90	102.499496	102.603837	102.616859	102.619028
180	102.499599	102.605953	102.617976	102.619681
360	102.500454	102.606351	102.618421	102.620069

Table 5. Table 5: Parameters of the model for the JP Morgan bond.

Parameters of the defaultable stock price model

$a_1=0.0312763$

$a_2=0.0356952$

$b_1=0.00038362$

$b_2=0.00172115$

$c=0.346622$

$\beta =-0.223027$

Parameters of the interest rate model

$\kappa =0.14485883018483803$

$\delta =0.01330207057173363$

$\theta =0.03467342840511061$

Correlation coefficient

$\rho =0.497108$

Initial conditions

$S_0=1.0$

$r_0=0.01469383913023823$

Table 6. Table 6: Value of the JP Morgan bond for different meshes and time steps.

Time steps per year	Mesh 4	Mesh 8	Mesh 16	Mesh 32
90	103.725041	103.596891	103.572191	103.570225
180	103.841153	103.605155	103.575270	103.572747
360	103.794567	103.602389	103.576483	103.574147

Table 7. Table 7: Comparison of PDE numerical solution with 360 time steps and Mesh 32 with respect to Monte Carlo 95 % percent 95 95\% confidence interval with 100000 simulations and 2nd-order asymptotic approximation. Note that the real market prices of the bonds of the UBS and JPM are respectively 102.62 and 103.57.

	PDE solution	MC confidence interval	asymptotic approximation
UBS	102.620069	[102.52477413, 102.72696887]	102.566
JPM	103.574147	[103.55570424, 103.65668368]	103.584

Equations

$$d{r}_t=\kappa (\theta -r_t)dt+\delta d{W}_t^1,$$

$$\sigma (t,S_t)=a(t)S_t^{\beta },$$

$$\lambda (t,S_t)=b(t)+c\sigma (t,S_t{)}^2=b(t)+ca(t{)}^2{S}_t^{2\beta },$$

$$d{S}_t=(r_t+\lambda (t,S_t))S_{t}dt+\sigma (t,S_t)S_{t}d{W}_t^2.$$

$$\tilde{\xi }=\text{inf}\left\{t\ge 0:{\int }_0^t\lambda (t,S_t)\ge e\right\},$$

$$V(0,S_0,r_0;T)=FV[\sum _{i=1}^{M}c{p}_i\mathbb{E}\left[\exp \left(-{\int }_0^{t_i}\left(r_u+\lambda (u,S_u)\right)du\right)\right]+\mathbb{E}\left[\exp \left(-{\int }_0^T\left(r_u+\lambda (u,S_u)\right)du\right)\right]$$

$$+\eta \left(1-\mathbb{E}\left[\exp \left(-{\int }_0^T\left(r_u+\lambda (u,S_u)\right)du\right)\right]-{\int }_0^T\mathbb{E}\left[\exp \left(-{\int }_0^{{\tau }_1}\left(r_u+\lambda (u,S_u)\right)du\right)r_{{\tau }_1}\right]d{\tau }_1\right)],$$

$$u_1(0,S_0,r_0;t_i)$$

$$u_1(0,S_0,r_0;T)$$

$$u_2(0,S_0,r_0;{\tau }_1)$$

$$V(0,S_0,r_0;T)$$

$$\left\{\begin{array}[]{ll}\overline{\mathcal{L}}\left[u(t,S,y)\right]=0,&t<T_{1},\,(S,y)\in(0,\,\infty)\times(-\infty,\,\infty),\\ u(T_{1},S,y)=h(S,y),&(S,y)\in(0,\,\infty)\times(-\infty,\,\infty),\end{array}\right.$$

$$\overline{\mathcal{L}}[u]$$

$$\left\{\begin{array}[]{ll}\overline{S}_{t}=S_{0}\exp(X_{t}){\mathbbm{1}}_{\{\xi>t\}}\qquad S_{0}>0\\ dr_{t}=\kappa(\theta-r_{t})dt+\delta\d{W}_{t}^{1}\\ dX_{t}=(r_{t}-\frac{1}{2}\sigma(t,X_{t})^{2}+\lambda(t,X_{t}))dt+\sigma(t,X_{t})dW_{t}^{2}\\ dW^{1}_{t}dW^{2}_{t}=\rho dt,\end{array}\right.$$

$$\mathcal{L}$$

$$=$$

$$\Sigma \left(t,x,r\right)$$

$$\gamma \left(t,x\right)$$

$${\mathcal{O}}_R:=\left\{(x,r)\in {\mathbb{R}}^2\mid \left|(x,r)\right|<R\right\};$$

$$u(t,x,r):=\mathbb{E}\left[h(X_T,r_T)\right]={\int }_{{\mathcal{O}}_R}\Gamma (t,x,r;T,z,s)h(z,s)dzds$$

$$\{\begin{array}{ll}\mathcal{L}u(t,x,r)=0,& (t,x,r)\in [0,T)\times {\mathcal{O}}_R,\\ u(T,x,r)=h(x,r)& (x,r)\in {\mathcal{O}}_R.\end{array}$$

$$x_0=t,{\tilde{x}}_1=S\text{and}{\tilde{x}}_2=y.$$

$${\Omega }^{\ast }=(0,x_0^{\infty })\times \left(\frac{1}{{\tilde{x}}_1^{\infty }},{\tilde{x}}_1^{\infty }\right)\times (-{\tilde{x}}_2^{\infty },{\tilde{x}}_2^{\infty })$$

$$x_1={\tilde{x}}_1-\frac{1}{{\tilde{x}}_1^{\infty }},x_2={\tilde{x}}_2+{\tilde{x}}_2^{\infty },$$

$$\tilde{\Omega }=(0,x_0^{\infty })\times (0,x_1^{\infty })\times (0,x_2^{\infty })$$

$${\tilde{\Gamma }}_i^{-}=\{(x_0,x_1,x_2)\in \tilde{\Gamma }\mid x_i=0\},{\tilde{\Gamma }}_i^{+}=\{(x_0,x_1,x_2)\in \tilde{\Gamma }\mid x_i=x_i^{\infty }\},i=0,1,2.$$

$$\mathcal{L}[u]=\sum _{i,j=0}^2{b}_{ij}\frac{{\partial }^{2}u}{\partial x_i{x}_j}+\sum _{j=0}^2{b}_j\frac{\partial u}{\partial x_j}+b_{0}u,$$

$$\mathbf{B}$$

$$\mathbf{b}$$

$$b_0$$

$${\Sigma }^0=\left\{x\in \tilde{\Gamma }/\sum _{i,j=0}^2{b}_{ij}{m}_i{m}_j=0\right\},{\Sigma }^1=\tilde{\Gamma }-{\Sigma }^0,$$

$${\Sigma }^2=\left\{x\in {\Sigma }^0/\sum _{i=0}^2\left(b_i-\sum _{j=0}^2\frac{\partial b_{ij}}{\partial x_j}\right)m_i<0\right\}.$$

$${\Gamma }_i^{-}=\{(x_1,x_2)\in \Gamma \mid x_i=0\},{\Gamma }_i^{+}=\{(x_1,x_2)\in \Gamma \mid x_i=x_i^{\infty }\},i=1,2.$$

$${\partial }_{\tau }u-Div(\mathbf{A}\nabla u)+\mathbf{v}\cdot \nabla u+lu=0$$

Full text

1 Introduction

Default risk can be modelled by using two different approaches: the structural approach and the reduced-form approach. For more information about the history, advantages and drawbacks of each one of the two approaches we refer the readers to [1, 8, 9]. In the structural method default occurs when the firm asset value reaches some lower barrier, whereas in the reduced-form approach, the default event is assumed to be unpredictable and governed by a default intensity process that could be either deterministic or stochastic. Moreover, in the reduced-form approach the default event can occur without any correlation with the firm value. In this work we consider the reduced-form approach.

In the literature there are several papers devoted to default risk modelling by using the two different approaches. On one hand, the first approach is taking into consideration in [1] where the author proves an exact formula for the valuation of defaultable coupon bonds by generalizing the one derived in [14]. On the other hand, both approaches are combined in [9] also to price defaultable bonds.

The main objective of this paper is to obtain the price of a defaultable coupon bond under the extended Jump to Default Constant Elasticity of Variance (JDCEV) model proposed in [7], also considering the possibility of incorporating negative interest rates which is introduced in [11]. The JDCEV model was first introduced in [7] as a hybrid credit and equity model, although only taking into account constant positive interest rates. In this work we assume that the interest rates dynamics is governed by a stochastic process which takes into account the possibility of negative interest rates as in [11]. The incorporation of negative interest rates results more realistic according to the current situation of real markets. As it was pointed out before, this model can be set in the framework of the reduced-form approaches. More precisely, in the present paper we define the instantaneous volatility of the stock price as a constant elasticity of variance (CEV) process and we assume that the default intensity is an affine function of the instantaneous variance of the underlying stock.

From the mathematical point of view, the valuation problem of a defaultable coupon bond can be posed in terms of a sequence of partial differential equation (PDE) problems, where the underlying stochastic factors are the interest rates and the stock price. Moreover, the stock price follows a diffusion process interrupted by a possible jump to zero (default), as it is indicated in [7]. In order to compute the value of the bond we need to solve two partial differential equation problems for each coupon, with maturities equal to those coupon payment dates. Concerning the numerical solution of these PDE problems, after a localization procedure to formulate the problems in a bounded domain and the study of the boundaries where boundary conditions are required following the ideas introduced in [17], we propose appropriate numerical schemes based on a Crank-Nicolson semi-Lagrangian method for time discretization combined with biquadratic Lagrange finite elements for space discretization. The numerical analysis of this Lagrange-Galerkin method has been addressed in [3, 4]. Once the numerical solution of the PDEs is obtained, a kind of post-processing is carried out in order to achieve the value of the bond. This post-processing includes the computation of an integral term which is approximated by using the composite trapezoidal rule.

This paper is organized as follows. In Section 2, we describe the two stochastic factors (i.e. the interest rate and the defaultable stock price) that are involved in the model, and we state the PDE problem that governs the valuation of non callable defaultable coupon bonds. In Section 3 we establish the existence of solution for the PDE problems. In Section 4, we formulate the pricing problem in a bounded domain after a localization procedure and we impose appropriate boundary conditions. Then, we introduce the discretization in time of the problem by using a Crank-Nicolson characteristic scheme, and we state the variational formulation of the problem in order to apply finite elements for the discretization in the asset and interest rate variables. In Section 5 we present some numerical results to illustrate the good performance of the models and numerical methods, also including a comparison with the results obtained with an alternative Monte Carlo technique and an asymptotic expansion method. Finally, we finish with some conclusions in Section 6.

2 Mathematical modelling

2.1 Stochastic underlying variables

A non callable defaultable coupon bond is a financial derivative product, the underlying variables of which are the interest rate and the defaultable stock price. The interest rate at time $t$, $r_{t}$, is assumed to be stochastic and its dynamics under a risk neutral probability measure is driven by the Vasicek model [20], in which the spot rate is governed by the Ornstein-Uhlenbeck process:

[TABLE]

where $\kappa>0$ is the speed of adjustment in the mean reverting process, $\theta>0$ is the long-term mean of the short-term interest rate, $\delta>0$ is the interest rate volatility and $W_{t}^{1}$ is the standarized Wiener process for the interest rate. Negative interest rates are taken into consideration with this model contrary to other interest rate dynamics, such as CIR model [10]. Note that the possibility of negative interest rates results more realistic according to the current situation of the markets.

In the JDCEV model we assume that the local volatility of the stock is given by

[TABLE]

where $S_{t}$ denotes the stock price at time $t$, $\beta<0$ is the elasticity parameter and $a(t)>0$ is the time dependent volatility scale function. Moreover, the default intensity can be written in terms of the stock volatility and the stock price in the following way

[TABLE]

where $b(t)\geq 0$ is a deterministic non-negative function of time and $c>0$ governs the sensitivity of the default intensity with respect to the volatility, as it is indicated in [7]. The other source of uncertainty, the predefaultable stock price at time t, $S_{t}$, under a risk neutral probability measure satisfies the following stochastic differential equation:

[TABLE]

where $W_{t}^{2}$ is the standarized Wiener process for the stock price. Both Wiener processes, $W_{t}^{1}$ and $W_{t}^{2}$ can be correlated according to $dW_{t}^{1}dW_{t}^{2}=\rho dt$, where $\rho$ denotes the instantaneous correlation coefficient such that $\mid\rho\mid<1$. Moreover, in this model the time of default $\xi$ has two parts. On one hand, we can define the predictable part as $\xi_{0}=\text{inf}\{t\geq 0:S_{t}=0\}$, whereas on the other hand the random part $\tilde{\xi}$ is given by the following expresion:

[TABLE]

where $e$ is an exponential random variable $e\sim Exp(1)$. Thus, the time of default is defined as $\xi=\xi_{0}\wedge\tilde{\xi}$.

Finally, the defautable stock price is given by $\overline{S}_{t}=S_{t}{\mathbbm{1}}_{\{\xi>t\}}$.

2.2 PDE formulation

In this work we consider the valuation of a non callable defaultable bond with the possibility of paying a series of coupons at given dates $t_{i}$, for $i=1,...,M$, where $M$ is the number of coupons and $T=t_{M}$ is the maturity of the bond. At maturity, the bond holder receives the face value (FV) plus the last coupon. Thus, let us denote by $cp_{i}$ the amount of money paid at coupon payment date $t_{i}$. This amount is computed by multiplying the coupon rate, the frequency of coupon payments and the FV. Having this in view, the value of a non callable defaultable coupon bearing bond at time $t=t_{0}=0$ for the spot values $S_{0}$ and $r_{0}$, $V(0,S_{0},r_{0};T)$, which is understood as the discounted value of the future coupon payments and the face value of the bond, is given by

[TABLE]

where $\eta$ is the recovery rate in case of default. An analysis of the recovery rate by industrial sector and debt seniority is carried out in [2].

Next, if we denote by

[TABLE]

then the expression of the bond value (2.2) can be written equivalently as

[TABLE]

Moreover, by applying the Feynman-Kac formula (see [19], for example) and using the change of variable $y_{t}=r_{t}\exp(\kappa t)$, the functions $u_{1}$ and $u_{2}$ are the solutions of the Cauchy problem

[TABLE]

with $h(S,y)=1$ for $u=u_{1}$ or $h(S,y)=\exp(-\kappa T_{1})y$ for $u=u_{2}$, respectively. Moreover, the operator $\overline{\mathcal{L}}$ is defined as follows

[TABLE]

3 Existence of solution

In this section, we aim to prove the existence of a “local” solution to the Cauchy problem (7). For this purpose, let $X$ be the logarithm of the pre-defaultable stock price, i.e. $X_{t}=\log{S_{t}}$, $t\in[0,T]$. Our model becomes

[TABLE]

with $\sigma(t,x)=a(t)e^{\beta x}$ and $\lambda(t,x)=b(t)+c\sigma(t,x)^{2}$. The corresponding infinitesimal generator is the operator $\mathcal{L}$ defined as

[TABLE]

where

[TABLE]

Operator $\mathcal{L}$ is only locally uniformly parabolic in the sense that, for any ball

[TABLE]

the coefficients of $\mathcal{L}$ satisfy the following conditions:

(H1) The matrix $\Sigma(t,x,r)$ is positive definite, uniformly with respect to $(t,x,r)\in\left(0,T\right]\times\mathcal{O}_{R}$.

(H2) The coefficients $\Sigma,\mu,\gamma$ are bounded and Hölder-continuous on $\left(0,T\right]\times\mathcal{O}_{R}$.

Under these conditions we can resort on the recent results in [18], Theor. 2.6, or [16], Theor.1.5, about the existence of a local density for the process $(X,r)$.

Theorem 3.1

For any $R>0$, the process $(X,r)$ has a local transition density on $\mathcal{O}_{R}$, that is a non-negative measurable function $\Gamma=\Gamma(t,x,r;T,z,s)$ defined for any $0<t<T$, $(x,r)\in{\mathbb{R}}^{2}$ and $(z,s)\in\mathcal{O}_{R}$ such that, for any continuous function $h=h(x,r)$ with compact support in $\mathcal{O}_{R}$, we have

[TABLE]

and $u$ satisfies

[TABLE]

Problem (18) can be used for numerical approximation purposes. However, notice that (18) does not have a unique solution due to the lack of lateral boundary conditions. Nevertheless, numerical schemes can be implemented imposing artificial boundary conditions and the result which guarantees the validity of such approximations is the so-called principle of not feeling the boundary. A rigorous statement of this result can be found in [13], Appendix A, or in [18], Lemma 4.11.

4 Numerical methods

In order to obtain a numerical approach of the value of a non callable defaultable coupon bond we need to solve the Cauchy problem (7) for $u=u_{1}$ and $u=u_{2}$ with maturity $T_{1}=t_{i}$ for $i=1,...,M$, that is each coupon payment date for both cases. Once these problems are solved, the value of the bond is given by expression (6) in which an integral term appears. This integral term will be approximated by means of the classical composite trapezoidal rule. For the numerical solution of the PDE problem, we propose a Crank-Nicolson characteristics time discretization scheme combined with a piecewise bi-quadratic Lagrange finite element method. The convergence properties of this Lagrange-Galerkin method have been mathematically analyzed in [3, 4] for time and space discretization. More recently, it has been applied to the valuation of pension plans without and with early retirement in [5] and [6], respectively. In order to apply this set of numerical techniques, first a localization procedure is used to cope with the initial formulation in an unbounded domain.

4.1 Localization procedure and formulation in a bounded domain

In this section we replace the unbounded domain by a bounded one. In order to determine the required boundary conditions for the associated PDE problem we follow [17] which is based on the theory proposed by Fichera in [12]. Let us introduce the following notation:

[TABLE]

For this purpose, let us consider both $\tilde{x}_{1}^{\infty}$ and $\tilde{x}_{2}^{\infty}$ to be large enough suitably chosen real numbers such that the solution in the region of financial interest is not affected by the truncation of the domain. Let

[TABLE]

with $x_{0}^{\infty}=T_{1}$. Additionally, we make the changes of variables

[TABLE]

as well as the values $x_{1}^{\infty}=\tilde{x}_{1}^{\infty}-\frac{1}{\tilde{x}_{1}^{\infty}}$ and $x_{2}^{\infty}=2\tilde{x}_{2}^{\infty}$, which leads to replace the bounded domain $\Omega^{*}$ by the following one:

[TABLE]

Then, let us denote the Lipschitz boundary by $\tilde{\Gamma}=\partial\tilde{\Omega}$ such that $\tilde{\Gamma}=\bigcup_{i=0}^{2}(\tilde{\Gamma}_{i}^{-}\cup\tilde{\Gamma}_{i}^{+})$, where

[TABLE]

Then, the operator defined in (8) can be written in the form:

[TABLE]

where the involved data are given by

[TABLE]

Thus, following [17], in terms of the normal vector to the boundary pointing inward $\tilde{\Omega}$, $\mathbf{m}=(m_{0},m_{1},m_{2})$, we introduce the following subsets of $\tilde{\Gamma}$:

[TABLE]

[TABLE]

As indicated in [17] the boundary conditions at $\Sigma^{1}\bigcup\Sigma^{2}$ for the so-called first boundary value problem associated with (20) are required. Note that $\Sigma^{1}=\{\tilde{\Gamma}_{1}^{-},\,\tilde{\Gamma}_{1}^{+},\,\tilde{\Gamma}_{2}^{-},\,\tilde{\Gamma}_{2}^{+}\}$ and $\Sigma^{2}=\{\tilde{\Gamma}_{0}^{+}\}$. Therefore, in addition to a final condition (see section 2.2), we need to impose boundary conditions on $\tilde{\Gamma}_{1}^{-}$, $\tilde{\Gamma}_{1}^{+}$, $\tilde{\Gamma}_{2}^{-}$ and $\tilde{\Gamma}_{2}^{+}$. Next, we will impose Dirichlet conditions on $\tilde{\Gamma}_{1}^{-}$, $\tilde{\Gamma}_{2}^{-}$ and $\tilde{\Gamma}_{2}^{+}$, whereas on $\tilde{\Gamma}_{1}^{+}$ we will impose a homogeneous Neumann condition.

Taking into account the previous change of spatial variable and making the change of time variable $\tau=T_{1}-t$, we write the equation (7) in divergence form in the bounded spatial domain $\Omega=(0,x_{1}^{\infty})\times(0,x_{2}^{\infty})$. Moreover, we decompose the boundary $\partial\Omega=\bigcup_{i=1}^{2}({\Gamma}_{i}^{-}\cup\Gamma_{i}^{+})$ of the spatial domain $\Omega$ as follows:

[TABLE]

Thus, the initial boundary value problem (IBVP) takes the following form:

Find $u:[0,T_{1}]\times\Omega\rightarrow\mathbb{R}$ such that

[TABLE]

For problem (27)-(30), the diffusion matrix $\mathbf{A}$ and the velocity field $\mathbf{v}$ are given by

[TABLE]

while the reaction function $l$, the initial condition $g$ and the function $f$ for the Dirichlet boundary conditions are defined as follows:

[TABLE]

4.2 Time discretization

The method of characteristics is based on a finite differences scheme for the discretization of the material derivative, i.e., the time derivative along the characteristic lines of the convective part of the equation (27). The material derivative operator is given by

[TABLE]

For a brief description of the method, we first define the characteristics curve through $\mathbf{x}=(x_{1},x_{2})$ at time $\bar{\tau}$, $X(\mathbf{x},\bar{\tau};s)$, which satisfies:

[TABLE]

In order to discretize in time the material derivative in equation (27), let us consider a number of time steps $N$, the time step $\Delta\tau=T/N$ and the time mesh points $\tau^{n}=n\Delta\tau,$ $n=0,\frac{1}{2},1,\frac{3}{2},\ldots,N$.

The material derivative approximation by the characteristics method for both problems is given by:

[TABLE]

where $u=u_{1},u_{2}$ and $X^{n}(\mathbf{x})=X(\mathbf{x},\tau^{n+1};\tau^{n})$. In this case, the solution of (33) is not computed analytically. Instead, we consider numerical ODE solvers to approximate the characteristics curves (see [3], for example). More precisely, in this work we employ the explicit second order Runge-Kutta method to approximate the values of $X^{n}(\mathbf{x})$.

Next, we consider a Crank-Nicolson scheme around $\left(X(\textbf{x},\tau^{n+1};\tau),\tau\right)$ for $\tau=\tau^{n+\frac{1}{2}}$. So, the time discretized equation for $u=u_{1},u_{2}$ can be written as follows:

Find $u^{n+1}$ such that:

[TABLE]

In order to obtain the variational formulation of the semi-discretized problem, we multiply (4.2) by a suitable test function, integrate in $\Omega$, use the classical Green formula and the non classical following one [15]:

[TABLE]

Note that, as the characteristics curves cannot be obtained analytically, the terms $(\nabla X^{n})^{-1}(\mathbf{x})$ and $Div((\nabla X^{n})^{-T}(\mathbf{x}))$ in (4.2) are replaced by the following approximations (see [3] for more details):

[TABLE]

[TABLE]

where $\mathbf{L}=\nabla\mathbf{v}$.

After the previous steps, we can write a variational formulation for the time discretized problem as follows:

Find $u^{n+1}\in H^{1}(\Omega)$ satisfying the Dirichlet boundary condition (21), such that:

[TABLE]

for all $\psi\in H^{1}(\Omega)$ such that $\psi=0$ on $\Gamma_{1}^{-}$, $\Gamma_{2}^{-}$ and $\Gamma_{2}^{+}$, In the last term $a_{12}$ is the corresponding coefficient of the diffusion matrix $\mathbf{A}$.

4.3 Finite elements discretization

For the spatial discretization we consider $\{\tau_{h}\}$, a quadrangular mesh of the domain $\Omega$. Let $(T_{2},\mathcal{Q}_{2},\Sigma_{T_{2}})$ be a family of piecewise quadratic Lagrangian finite elements, where $\mathcal{Q}_{2}$ denotes the space of polynomials defined in $T_{2}\in\tau_{h}$ with degree less or equal than two in each spatial variable and $\Sigma_{T_{2}}$ the subset of nodes of the element $T_{2}$. More precisely, let us define the finite elements space $u_{h}$ by

[TABLE]

where $\mathcal{C}^{0}(\Omega)$ is the space of piecewise continuous functions on $\Omega$.

4.4 Composite trapezoidal rule

In order to obtain the value of the bond at origination, i.e. $V(0,S_{0},r_{0};T)$, by means of expression (6) the computation of an integral term is required. The approximation of this integral is carried out by using a suitable numerical integration procedure. More precisely, we employ the classical composite trapezoidal rule with $M+1$ points, where $M$ is the number of coupons, in the following way:

[TABLE]

where $h=\frac{T}{M}$, $k_{j}=jh$ for $j=1,...,M-1$ and $u_{2}(0,S_{0},r_{0};0)=r_{0}$.

5 Numerical results

In order to show the good performance of the model and numerical methods explained in Section 4, we present some numerical results. In the following examples, the value of some of the parameters involved in the underlying factors are taken from the literature, in particular from [11] where the authors calibrate the model parameters to market data. More precisely, first the interest rate model is calibrated to zero-coupon bonds (ZCB) and next the model is calibrated to CDS spreads, the price of which is obtained by means of an asymptotic expansion method.

In both examples, the number of elements and nodes of the finite element meshes employed in the numerical solution of the problems are shown in Table 1.

5.1 Example 1

First, we consider the simple case of the valuation of default-free zero-coupon bonds with different maturities. In this setting, the valuation problem is reduced to a one-factor model. The purpose of this example is to compare the value of the bonds we obtain with the market zero-coupon curve. In order to obtain the value of the bonds, we solve the IBVP (27)-(22) with initial condition $g(x_{1},x_{2})=1$ (corresponding to $h(S,r)=1$)and only taking into account as underlying factor the interest rate. The value of the parameters involved in the interest rate model are the ones collected in Table 2. The values of the zero-coupon curve and the approximated ones obtained by solving the here proposed model are presented in Table 3. For computing the numerical solution of Table 3 we consider Mesh 32 from Table 1 and the time step $\Delta\tau=\frac{1}{360}$ (one day).

5.2 Example 2

Next, we consider the valuation of two real defaultable bonds traded in the market and issued by different firms. For this purpose, as in [11], we assume that the coefficients $a(t)$ and $b(t)$ in (2) are linearly dependent on time and their expresions are given by $a(t)=a_{1}t+a_{2}$, $b(t)=b_{1}t+b_{2}$, where $a_{1}$, $a_{2}$, $b_{1}$ and $b_{2}$ are constants.

On one hand, we take into account the pricing of a bond from UBS with maturity 5 years and a face value of 100. The bond pays annually coupon rates of 1.25 basis points and the recovery rate at the event of default is 40$\%$. The model parameter values for this example are collected in Table 2. In this case the correlation coefficient $\rho$ is assumed to be zero. Next, in Table 4 we present the value of the bond for different meshes and time steps. In this case, we can appreciate that the price of the bond converges to 102.62.

On the other hand, we present a correlated case. More precisely, we address the valuation of a bond from JP Morgan with maturity 5 years and a face value of 100. The bond pays annually coupon rates of 3.25 basis points and the recovery rate at the event of default is again 40$\%$. For this second example the parameter values of the model are the ones which appear in Table 5. As we have just pointed out the correlation coefficient $\rho$ is different from zero. Finally, the value of this bond is shown in Table 6. In this case, the value of the bond converges to 103.57.

In both cases, in order to obtain the value of the bond for we need to solve for each coupon payment date the IBVP (27)-(22) with maturity equal to those dates and with initial condition $g(x_{1},x_{2})=1$ (corresponding to $h(S,r)=1$) or $g(x_{1},x_{2})=\exp(-\kappa T)(x_{2}-\tilde{x}_{2}^{\infty})$ (corresponding to $h(S,r)=r$. More precisely, in both examples the maturity of the bond is 5 years and the frequency of coupon payments is annually, thus to obtain the value of the bond we need to solve the problem (27)-(22) 10 times, i.e. 5 times with one initial condition to obtain the value of $u_{1}$ and 5 times with the other initial condition to obtain the value of $u_{2}$.

Next, in Figures 1 and 2 we show the mesh value of the UBS bond and the JP Morgan bond, respectively. Both figures are obtained with the finest mesh and with time step $\Delta\tau=\frac{1}{360}$ (one day).

Finally, for UBS and JPM bonds the Table 7 shows a comparison between the results obtained with the proposed numerical method for the PDE model, a crude Monte Carlo technique and the asymptotic approximation method introduced in [21, 22]. More precisely, for the PDE numerical results we consider Mesh 32 with $360$ time steps while we use $100000$ simulations in Monte Carlo and we show the $95\%$ confidence interval.

In the fourth column, we present the prices of the UBS and JPM bonds computed with the asymptotic approximation method. This method consists of approximating the solution $u$ of the parabolic PDE (18) by applying Theorem 3.2 in [21]. Note that the obtained value from the PDE numerical methods belongs to the confidence interval and is more accurate than the asymptotic approximation when compared with the real market prices.

Concerning computational time, the asymptotic method takes around 0.047 seconds to obtain one value and Monte Carlo simulation takes 74 seconds for UBS bond and 45 seconds for JPM bond to obtain one value, while the numerical solution of the PDE takes around 440 seconds to obtain the bond values at all the 4225 mesh nodes.

6 Conclusions

In this paper we have considered the valuation of a non callable defaultable coupon bond where the underlying stochastic factors are the interest rate and the defaultable stock price. The pricing problem is posed as a sequence of IBVPs. More precisely, two PDE problems with different initial conditions with maturity each coupon payment date need to be solved. Once the numerical solution of these problems is carried out, the value of the bond is computed by means of an expression which also involves the computation of an integral term.

In order to obtain a numerical solution of the PDE problems, we have proposed appropriate numerical methods based on Lagrange-Galerkin formulations. More precisely, we combine a Crank-Nicolson semi-Lagrangian scheme for time discretization with biquadratic Lagrange finite elements for space discretization. Moreover, the integral term which is involved in the computation of the bond value is approximated by means of the classical composite trapezoidal rule. Finally, we show some numerical results in order to illustrate the behaviour of the proposed methods and its a agreement with the ones obtained by an alternative Monte Carlo technique and an asymptotic aproximation method.

7 Acknowledgements

First and fourth authors have been partially supported by the Spanish government (Ministerio de Economía y Competitividad, project MTM2016-76497-R) and Xunta de Galicia (grant GRC2014/044, including FEDER funds). Second, third and fourth authors have been supported by EU H2020-MSCA-ITN-2014 (WAKEUPCALL Grant Agreement 643045).