An iterative method for solving a PDE with free boundary arising from pricing corporate bond with credit rating migration

1.   Introduction

In recent years, with the frequent occurrence of financial risk events, more and more attention has been paid to the credit risks of financial products. Credit risks of financial products include both default risks and credit grade migration risks. The previous research pays more attention to default risk, but now credit grade migration risk has become an important role in the bond risk managements. The upgrade or downgrade of credit rating will affect the value of corporate bond. The free boundary models have been established in ^{[2,6,8,10,11,12,15]} for pricing corporate bonds with the characteristic of credit grade migration risk, in which the free boundary is determined by the ratio of corporate debt to corporate value.

In this paper we study the following PDE with free boundary for pricing a corporate bond with the characteristic of credit grade migration risk ^[6,13]

with the final value condition at the expiration time $T$

Here $S$ is the corporate asset value, $\tau$ is the time, $v_{L}(S, \tau)$ and $v_{H}(S, \tau)$ are the bond values in low and high credit grades respectively, $\sigma_{L}$ and $\sigma_{H} \ (0 < \sigma_{H} < \sigma_{L})$ are volatilities of the corporate asset value under the low and high credit rating respectively, $\gamma \ (0 < \gamma < 1)$ is the threshold ratio of corporate debt to corporate asset value, $r$ is the risk-free rate of interest, and $F$ is the face value of the bond. Generally, it can be assumed that $F = 1$. The defined domain of the free boundary problem is divided into a low rating region $\Omega_{L}$ where $0 < S < \frac{1}{\gamma}v_{L}$ and a high rating region $\Omega_{H}$ where $S > \frac{1}{\gamma}v_{H}$. It has been proved that two domains are separated by a free boundary $s(\tau)$, and

At the credit rating migration boundary $s(\tau)$, the values of the bond in low and high credit rating satisfy

where $s(\tau)$ is an apriorily unknown function since the solutions $v_{L}$ and $v_{H}$ are two apriorily unknown functions. It has been proved in [6,Theorems 5.1 and 6.1] that the free boundary problems (1.1)–(1.5) has a unique weak solution $\left(v(S, t), s(t)\right)$ with $v(S, t)\in W_{\infty}^{2, 1}\left(\left(\left(-\infty, \infty\right)\times [0, T]\right)\setminus \bar{Q}_{\rho}\right)\cap W_{\infty}^{1, 0}\left(\left(-\infty, \infty\right)\times [0, T]\right)$ and $s(t)\in C[0, T]$ for any $\rho > 0$, where

and $Q_{\rho} = \left(-\rho, \rho\right)\times \left(0, \rho^{2}\right)$.

The above problem requires not only solving the value of the bond, but also solving the free boundary. In financial engineering, pricing financial products with free boundary has long been recognized as a very challenging problem. A few numerical methods have been used to solve such problems. Explicit difference schemes are used in ^[6,11,15] to solve free boundary problems for pricing corporate bonds with credit grade migration risks. A front fixing method is derived in ^[7] to solve problems (1.1)–(1.5), which transforms the free boundary into a fixed boundary by including the unknown boundary into the equation, resulting in the differential equation becoming a nonlinear equation. For the transformed fixed boundary problem, the predictor-corrector algorithm and Newton-like iterative algorithms are used to solve the difference equations in ^[7]. The predictor-corrector algorithm is also an explicit discrete scheme that needs to satisfy the stability conditions, while the Newton-like iterative method needs a lot of computation to solve the nonlinear difference equations.

In this paper, we propose a novel method to solve the PDE with free boundary (1.1)–(1.5). An iterative algorithm is designed to generate weak solution sequences of fixed internal boundary problems. It is proved that both weak solution sequences are convergent. Since it is not easy to obtain analytical solutions of the fixed internal boundary problems, numerical methods are used to solve them. In each iteration step, an implicit-upwind difference scheme is applied to solve the fixed internal boundary problem. The stability and convergence order of the discrete scheme are given. Numerical experiments verify that the limit of the weak solution sequence is the solution of the free boundary problem and also verify that the discrete scheme is stable and first-order convergent. The advantages of this method are reflected in three aspects: first, the free boundary problem is transformed into a sequence of fixed internal boundary problems, which simplifies the problem and deepens the understanding of this free boundary problem; second, the implicit scheme is used to solve the fixed internal boundary problem in each iteration step, so as to ensure the stability of the discrete scheme without additional constraints; third, this method only involves solving the root of a single nonlinear equation without solving the system of nonlinear equations, which reduces the amount of calculation.

2.   Iterative method

By using the variable transformations $x = \ln S$ and $t = T-\tau$, and defining

we can derive the following equation from (1.1)–(1.5) and the assumption $F = 1$

where $x^{*}(t)$ is the free boundary transformed from $s(\tau)$. Here $x^{*}(t)$ is an apriorily unknown function since it should be solved by the following equation

where the solution $u$ is also an apriorily unknown function.

Let $u_{H}(x, t)$ and $u_{L}(x, t)$ be the solutions of problems

respectively. Then, using the method for solving the classical Black-Scholes equation ^[9], we can get the solutions of problems (2.3) and (2.4) as follows

where

and

Similar results have been given in the literature ^[12]. Furthermore, the following result can be obtained, which also has been proved in [12,Theorem 2.2].

Lemma 2.1 Let $u(x, t), u_{H}(x, t)$ and $u_{L}(x, t)$ be the solutions of problems (2.1), (2.3) and (2.4) respectively. Then, we have

Next, according to the theory of the linear parabolic equation ^[4], we construct two weak solution sequences $\left\{\bar{u}^{(k)}\right\}$ and $\left\{\underline{u}^{(k)}\right\}$ with $\bar{u}^{(k)}, \underline{u}^{(k)}\in W_{\infty}^{2, 1}\left(\left(\left(-\infty, \infty\right)\times [0, T]\right)\setminus \bar{Q}_{\rho}\right)\cap W_{\infty}^{1, 0}\left(\left(-\infty, \right.\right.$ $\left.\left.\infty\right)\times [0, T]\right)$ for $\rho > 0$ and $Q_{\rho} = \left(-\rho, \rho\right)\times \left(0, \rho^{2}\right)$, which are generated as follows

and

When the solution $\bar{u}^{(k-1)}$ of the $k-1$ iteration is known, the existence and uniqueness of $\bar{x}^{(k)}$ can be derived from the results in Lemmas 2.2 and 2.3. When $\bar{x}^{(k)}$ is known, the iteration equation of the $k$-th order is a fixed internal boundary problem. For the fixed internal boundary problems in (2.5), they are parabolic equations with discontinuous coefficients as discussed in ^[3,5,14]. By using a construction method as used in ^[5] we can prove that the parabolic equation with discontinuous coefficients exists a solution. Let $\bar{u}_{1}^{(k)}(x, t)$ and $\bar{u}_{2}^{(k)}(x, t)$ be particular solutions of the following differential equations respectively

and

Consider the following function

where $\phi_{1}(x, t)$ and $\phi_{2}(x, t)$ are the solutions of the following parabolic problems respectively

and

By imposing the conditions

we can get $A(t)$ and $B(t)$. From this we conclude that the parabolic equation with discontinuous coefficients in (2.5) exists a solution. Similar results can be obtained for (2.6).

Next, we give some properties of iterative solutions.

Lemma 2.2 The weak solutions $\bar{u}^{(k)}$ and $\underline{u}^{(k)}$ of problems (2.5) and (2.6) satisfy

for $k\geq 0$ and

for $k\geq 1$.

Proof. Hu et al. ^[6] regard the free boundary problem (2.1) as a parabolic equation with discontinuous coefficients, and apply the maximum principle to prove in Lemmas 4.2 and 4.5 and Theorem 5.1 of ^[6] that the properties (2.7)–(2.10) hold true for the solution of the free boundary problem (2.1). For the fixed internal boundary problems in (2.5) and (2.6), they are also parabolic equations with discontinuous coefficients. The only difference between the two equations is that the coefficient $\sigma$ is different. As long as $\sigma = \sigma_{H}+\left(\sigma_{L}-\sigma_{H}\right)H\left(u-\gamma e^{x}\right)$ in ^[6] is replaced by $\sigma = \sigma_{H}+\left(\sigma_{L}-\sigma_{H}\right)H\left(\bar{x}^{(k)}(t)-\gamma e^{x}\right)$ or $\sigma = \sigma_{H}+\left(\sigma_{L}-\sigma_{H}\right)H\left(\underline{x}^{(k)}(t)-\gamma e^{x}\right)$, it can be proved by the same method that the results (2.7)-(2.10) for the fixed internal boundary problems in (2.5) and (2.6) also hold true for $0 < t\leq T$ and the inequalities (2.7)–(2.10) become equations for $t = 0$, where $H(\xi)$ is the Heaviside function and $u$ is the solution of the free boundary problem.

In order to simplify the expression, we introduce the following problems

for $j = 1, 2$, which are any two iterative equations in the iterative problems (2.5) and (2.6). Furthermore, let $s_{j}(t)$ with $j = 1, 2$ be the solutions of the following problems

respectively. By making the variable transformation $v_{j} = e^{x}w_{j}$, problems (2.11) and (2.12) can be reduced to

and

for $j = 1, 2$, respectively.

Applying Lemma 2.2 we can get

and

for $j = 1, 2$. Combining the variable transformation $v_{j}(x, t) = e^{x}w_{j}(x, t)$ and inequalities (2.15) and (2.16) we obtain

and

Lemma 2.3 For each $j$, the problem (2.12) have a unique solution $s_{j}(t)$. Then each iterative equation in problems (2.5) and (2.6) has a unique solution.

Proof. For each $j$, it is assumed that there exist two solutions $s_{j}^{1}(t)$ and $s_{j}^{2}(t)$ to the problem (2.12). Suppose there exists $t_{0}\in [0, T]$ such that

Since $s_{j}^{1}(0) = s_{j}^{2}(0) = 0$, we have $t_{0}\neq 0$. Moreover, from (2.14) and (2.17) we have

which is a contradiction. Hence, for each $j$ we have $s_{j}^{1}(t)\leq s_{j}^{2}(t)$. Similarly, for each $j$ we also can get $s_{j}^{2}(t)\leq s_{j}^{1}(t)$. Therefore, for each $j$ we have $s_{j}^{1}(t) = s_{j}^{2}(t)$, which implies that the problem (2.12) have a unique solution $s_{j}(t)$ for each $j$.

Furthermore, it is easy to prove that each iterative equation in problems (2.5) and (2.6) has a unique solution by using the maximum principle as given in ^[3,5,14].

Lemma 2.4 For the solutions $v_{j}(x, t)$ of problems (2.11), the following results hold true:

(ⅰ) if $y_{1}(t)\geq y_{2}(t)$ for $t\in [0, T]$, then $v_{1}(x, t)\leq v_{2}(x, t)$ for $(x, t)\in \left(-\infty, \infty\right)\times [0, T]$ and $s_{1}(t)\leq s_{2}(t)$ for $t\in [0, T]$;

(ⅱ) if $y_{1}(t)\leq y_{2}(t)$ for $t\in [0, T]$, then $v_{1}(x, t)\geq v_{2}(x, t)$ for $(x, t)\in \left(-\infty, \infty\right)\times [0, T]$ and $s_{1}(t)\geq s_{2}(t)$ for $t\in [0, T]$.

Proof. Set $z(x, t) = v_{2}(x, t)-v_{1}(x, t)$ and

Then, if $y_{1}(t)\geq y_{2}(t)$ for $t\in [0, T]$, from (2.11) and (2.19) we have

which implies

where we have used Lemma 2.2. Obviously, $z(x, 0) = 0$. Hence, it follows by the maximum principle in the sense of weak solution that

which implies

By making the variable transformations $v_{j} = e^{x}w_{j}$, problems (2.11) can be reduced to problem (2.13). Set $\psi(x, t) = w_{2}(x, t)-w_{1}(x, t)$ and

Then, if $y_{1}(t)\geq y_{2}(t)$ for $t\in [0, T]$, from (2.13) and (2.21) we have

which implies

where we have used (2.18). It is also obvious that $\psi(x, 0) = 0$. Hence, it follows by the maximum principle in the sense of weak solution that

which implies

Furthermore, we can prove the following inequality

holds true. If (2.24) does not hold for all $(x, t)\in \left(-\infty, \infty\right)\times (0, T]$, then there is a point $\left(x_{0}, t_{0}\right)\in \left(-\infty, \infty\right)\times (0, T]$ such that $\psi(x, t)$ reaches its minimum value at $\left(x_{0}, t_{0}\right)$, i.e.,

Then, the contradiction can be drawn from (2.22) and (2.25). Hence, we have

Based on the inequality (2.26) we prove the following inequality

holds true. Since $s_{1}(0) = s_{2}(0) = -\ln \gamma$, (2.27) holds true for $t = 0$. If (2.27) is not valid for all $t$, then there exists $t_{0}\in (0, T]$ such that

Then, from (2.14), (2.17) and (2.26) we have

which is a contradiction.

Combining inequalities (2.20) and (2.27), we conclude that the results in (ⅰ) hold true.

Using the same method used in proving (ⅰ), we also can prove that the results in (ⅱ) hold true.

The following theorems give the convergence of iterative sequences $\left\{\bar{u}^{(k)}\right\}$ and $\left\{\underline{u}^{(k)}\right\}$ in (2.5) and (2.6), which are the main results of this paper.

Theorem 2.5 The solution sequence $\left\{\bar{u}^{(k)}, \bar{x}^{(k)}\right\}$ of problem (2.5) satisfies

which imply that the weak solution sequence $\left\{\bar{u}^{(k)}\right\}$ and the internal boundary sequence $\left\{\bar{x}^{(k)}\right\}$ are convergent respectively.

Proof. From Lemma 2.1 we have

Combining (2.30) with Lemma 2.4 we can obtain

Then, by using Lemma 2.4 and (2.31) we can get

Set $\bar{v}^{(1)} = \bar{u}^{(1)}-u_{L}$ and

Then, from (2.4) and (2.5) we have

which implies

where we have used Lemma 2.2. It is also obvious that $\bar{v}^{(1)}(x, 0) = 0$. Hence, it follows by the maximum principle in the sense of weak solution that $\bar{v}^{(1)}\geq 0$, i.e.,

From Lemma 2.4 and (2.32) we have

Then, by using Lemma 2.4 and (2.34) we can get

Set $\bar{v}^{(2)} = u_{H}(x, t)-\bar{u}^{(2)}$ and

Then, from (2.3) and (2.5) we have

which implies

where we also have used Lemma 2.2. It is also obvious that $\bar{v}^{(2)}(x, 0) = 0$. Hence, it follows by the maximum principle in the sense of weak solution that $\bar{v}^{(2)}\geq 0$, i.e.,

From (2.35), (2.36) and Lemma 2.4 we have

Then, by using Lemma 2.4 and (2.37) we can get

Thus, from (2.38) and Lemma 2.4 we can get

Furthermore, by using Lemma 2.4 and (2.39) we have

Next we assume that the inequalities (2.28) and (2.29) hold true when the number of iterations is not greater than $2k$. Then we have

Thus, from (2.41) and Lemma 2.4 we have

Furthermore, by using Lemma 2.4 and (2.42) we can get

Hence, it can be seen from the induction that the inequalities (2.28) and (2.29) hold true for all $k$. Thus, by using Arzelà-Ascoli Theorem, we can prove that the monotone bounded sequences $\left\{\bar{x}^{(2k-1)}\right\}$ and $\left\{\bar{u}^{(2k-1)}\right\}$ are convergent respectively. Similarly, we can also prove that the monotone bounded sequences $\left\{\bar{x}^{(2k)}\right\}$ and $\left\{\bar{u}^{(2k)}\right\}$ are convergent respectively.

Remark. Although Theorem 2.5 does not prove that the limit of the solution sequence is the solution of Eq (2.1) in the classical sense, it can be considered to prove that the limit satisfies Eq (2.1) in the sense of distribution ^[1]. Considering that this paper focuses on numerical calculation, we use numerical experiments to verify that the limit is the solution of Eq (2.1).

Using the same method for proving Theorem 2.5 we also can obtain the following results.

Theorem 2.6 The solution sequence $\left\{\underline{u}^{(k)}, \underline{x}^{(k)}\right\}$ of problem (2.6) satisfies

which imply that the weak solution sequence $\left\{\underline{u}^{(k)}\right\}$ and the internal boundary sequence $\left\{\underline{x}^{(k)}\right\}$ are convergent respectively.

3.   Discretization

Since it is difficult to get the analytical solution for problems (2.5) and (2.6), we use an implicit-upwind difference scheme to solve them.

First, the spatial domain $(-\infty, \infty)$ is truncated into a finite domain $[x_{\min}, x_{\max}]$. The boundary conditions are chosen to be $u(x_{\min}, t) = u_{L}\left(x_{\min}, t\right)$ and $u(x_{\max}, t) = u_{H}\left(x_{\max}, t\right)$. Generally, the error caused by the truncation of the domain is negligible for the value of the bond. A uniform mesh $\Omega^{N\times K} = \Omega^{N}\times \Omega^{K}$ is utilized to discretize the definition domain $[x_{\min}, x_{\max}]\times [0, T]$, where

and

For the differential operator

an upwind difference scheme on $\Omega^{N}$ is utilized to approximate the spatial derivatives and an implicit Euler method on $\Omega^{K}$ is utilized to approximate the time derivative:

where $\sigma = \sigma_{L}$ or $\sigma = \sigma_{H}$, and

For the differential operator

an upwind difference scheme also is utilized to approximate the left and right derivatives:

It is easy to know that the discrete scheme satisfies the maximum principle, which can be derived from the fact that the matrix related to the discrete operator $\left\{L^{N, K}, l^{N, K}\right\}$ is an M-matrix. Then, we can conclude that the discrete scheme is unconditionally stable and is first-order convergent by the maximum principle.

Furthermore, the nonlinear equation

in problems (2.5) and (2.6) can be solved by Newton iteration method, where $u$ can be approximated by numerical solutions. In general, the solution of the nonlinear equation does not happen to be the mesh point of $\Omega^{N}$. We choose the closest mesh point as the approximate solution of the nonlinear equation.

4.   Numerical experiments

In this section we present some numerical results to indicate experimentally the efficiency and accuracy of our method. We consider the same example as given in ^[6,7].

Example The PDEs (1.1)–(1.5) with parameters $\sigma_{L} = 0.4, \sigma_{H} = 0.2, r = 0.05, F = 1, \gamma = 0.8, \; T = 5, x_{\min} = -\ln 5, x_{\max} = \ln 5$.

To numerically calculate Eqs (2.5) and (2.6), the stopping criterion of the iterative algorithm is chosen as

where $\bar{U}$ and $\underline{U}$ are the numerical solutions of Eqs (2.5) and (2.6) respectively.

The comparison between our numerical results and those of the explicit difference method given in ^[6] shows that they are very consistent, which are presented in Table 1.

Figure 1 displays the iterative solutions at $t = 0$ for the iterative equation (2.5), which shows that the even number of iterative solutions are above the numerical solution $U$ and the odd number of iterative solutions are below the numerical solution $U$, and the iterative solutions are closer to the solution $U$ with the increase of the number of iterations. Figure 2 displays the iterative solutions at $t = 0$ for the iterative equation (2.6), which shows that the odd number of iterative solutions are above the numerical solution $U$ and the even number of iterative solutions are below the numerical solution $U$, and the iterative solutions are also closer to the solution $U$ with the increase of the number of iterations. Figure 3 gives the numerical solution of the corporate bond with credit rating migration and Figure 4 gives the numerical solution of the free boundary caused by credit rating migration. It's easy to see from Figures 3 and 4 that the function of the bond value has been decomposed into two regions by a free boundary and the free boundary is decreasing as expected.

Since the example has no analytical solution, the double-layer mesh principle is utilized to calculate the error and the corresponding convergence rate. The error in the discrete maximum norm is denoted by

and the convergence rate is denoted by

The maximum errors, convergence rates and number of iterations in the calculation of (2.5) and (2.6) for Example are presented in Table 2, which show that the discrete scheme is stable and first-order convergent.

From the perspective of convergence order, our numerical method and the methods given in ^[6,7] are first-order convergent. Compared with the explicit difference methods given in ^[6,7], we use the implicit scheme to solve the fixed boundary problem in each iteration step, which ensures the stability of our discrete scheme without additional constraints. Compared with the front fixing method with Newton-like iterative algorithms given in ^[7], our method only involves solving the root of a single nonlinear equation without solving the system of nonlinear equations, which reduces the amount of calculation.

5.   Conclusions and discussion

An iterative method for a PDE with free boundary arising from pricing corporate bond with credit grade migration risk has been proposed. The key to the success of this method is that the constructed iterative algorithm produces two weak solution sequences of fixed internal boundary problems which are proved to be convergent. Since it is not easy to obtain analytical solutions of the fixed internal boundary problems, numerical methods are used to solve them. In each iteration step, an implicit-upwind difference scheme is used to solve the fixed internal boundary problem, which ensures the stability of the discrete scheme without additional constraints. Moreover, this method only involves solving the root of a single nonlinear equation without solving the system of nonlinear equations, which reduces the amount of calculation. It is shown that the scheme is stable and first-order convergent. Numerical experiments verify that the limit of the weak solution sequence is the solution of the free boundary problem, and numerical experiments also verify the stability and convergence order of the discrete scheme. The study of this paper broadens the method of solving free boundary problems. In future we extend this method to solve two dimensional models ^[11].

Acknowledgments

We would like to thank the anonymous reviewers for their valuable suggestions and comments for the improvement of this paper. The work was supported by Major Humanities and Social Sciences Projects in Colleges and Universities of Zhejiang Province (Grant No. 2018GH020).

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this article.