Multiscale Analysis on the Pricing of Intensity-Based Defaultable Bonds

Abstract

This paper studies the pricing of intensity-based defaultable bonds where the volatility of default intensity is assumed to be random and driven by two different factors varying on fast and slow time scales. Corrections to the constant intensity of default are obtained and then how these corrections influence the term structure of interest rate derivatives is shown. The results indicate that the fast scale correction produces a more significant impact on the bond price than the slow scale correction and the impact tends to increase as time to maturity increases.

1. Introduction

In finance, the payoff for a defaultable bond would be less than the promised amount when the risky asset of a firm may default. The simplest type of the defaultable bond can be modeled by defining the time of default as the first arrival time of a Poisson process with a constant mean arrival rate called intensity. However, it becomes a common sense nowadays that the default intensity should be treated as a stochastic process depending on the macroeconomic environment. Refer to Jarrow and Protter [1], Duffie and Singleton [2], and Bielecki and Rutkowski [3] for general reference on default intensity. See also Jarrow and Turnbull [4], Lando [5], Schonbucher [6] and Musiela and Rutkowski [7] for major mathematical developments in default modeling.

This paper considers a stochastic intensity model of which motivation is described as follows. There is a recent paper by Liang et al. [8] that studied the limitation of methods for pricing credit derivatives under a reduced form framework if the default intensity process follows the Vasicek model [9]. In fact, the intensity given by the Vasicek model could be negative, whereas it should not be the case. So, this work adopts the Cox-Ingersoll-Ross (CIR) model [10] for the intensity of default. On the other hand, there are studies by B. Kim and J.-H. Kim [11] and Papageorgiou and Sircar [12] presenting the pricing of defaultable derivatives under a diffusion model for the default intensity. These works are based upon multiscale stochastic volatility described by Fouque et al. [13], where volatility follows both fast and slow scale variations. In particular, the authors of [12] showed an empirical evidence of the existence of two different scales by the calibration of the model on corporate yield curves. However, these papers consider a Vasicek model for the interest rate which severely limits the practical applications of the results since the interest rate is always positive. So, we take both the underlying interest rate and the intensity of default as CIR models whose solutions are always positive. In this sense, our model is a fundamental extension of the aforementioned models.

Also, as shown in [14], modeling the intensity of default in terms of a Cox process turned out to be inappropriate for producing loss distributions to capture real data which exhibits a heavier tail. However, as stated by Papageorgiou and Sircar [15] on multiname credit derivatives, the introduction of multiscale stochastic volatility in the default intensity process is enough to allow for a heavier tail in the loss distribution, which is compatible with empirical evidence from real data and also offsets the need for jump characteristics and maintains closed-form expressions for the conditional loss distribution. Given this observation, the volatility of the CIR model for default intensity is assumed to be given by a function of two different time-scale stochastic factors. The main concern of this paper is to investigate how these factors influence the interest rate derivatives.

The remaining sections are organized as follows. In Section 2, dynamics of a defaultable bond are formulated in terms of a system of multiscale stochastic differential equations of the CIR type and transformed into an asymptotic partial differential equation. In Section 3, the solution of it is approximated by using asymptotic analysis and subsequently the convergence error is estimated. Section 4 studies, numerically, the impact of the two scale factors in the stochastic intensity of default on the price and the subsequent yield curve of the bond. Section 5 is devoted to the pricing of a European option for a defaultable bond. Final remarks are given in Section 6.

2. Problem Formulation

In terms of short rate 𝑟_𝑡, intensity 𝜆_𝑡, two small positive parameters 𝜖 and 𝛿, and two processes 𝑌_𝑡 and 𝑍_𝑡 representing a fast scale factor and a slow scale factor of the volatility of the intensity 𝜆_𝑡, respectively, we assume that dynamics of the joint process (𝑟_𝑡,𝜆_𝑡,𝑌_𝑡,𝑍_𝑡) are given by the following stochastic differential equations (SDEs): 𝑑𝑟_𝑡=𝑎𝑟^∗−𝑟_𝑡𝑑𝑡+𝜎√𝑟_𝑡𝑑𝑊^𝑟_𝑡,𝑑𝜆_𝑡=𝑎𝜆^∗−𝜆_𝑡𝑑𝑡+𝑓𝑌_𝑡,𝑍_𝑡𝜆_𝑡𝑑𝑊^𝜆_𝑡,𝑑𝑌_𝑡=1𝜖𝑚−𝑌_𝑡𝜈√𝑑𝑡+2√𝜖𝑑𝑊^𝑦_𝑡,𝑑𝑍_𝑡=𝛿𝑔𝑍_𝑡𝑑𝑡+√𝛿ℎ𝑍_𝑡𝑑𝑊^𝑧_𝑡 (1) under a risk-neutral probability measure (or equivalent martingale measure), where 𝑊^𝑟_𝑡, 𝑊^𝜆_𝑡, 𝑊^𝑦_𝑡, and 𝑊^𝑧_𝑡 are standard Brownian motions with a correlation structure given by the coefficients 𝜌_𝑟𝜆, 𝜌_𝑟𝑦, 𝜌_𝑟𝑧, 𝜌_𝜆𝑦, 𝜌_𝜆𝑧, and 𝜌_𝑦𝑧 with 𝜌_𝑟𝜆=0. The function 𝑓∶ℝ²→ℝ⁺ is assumed to be bounded, smooth, and strictly positive and the functions 𝑔 and ℎ satisfy the Lipschitz and growth conditions so that the corresponding SDE admits a unique strong solution. We note particularly that the process 𝑌_𝑡 is an ergodic process whose invariant distribution is given by the Gaussian probability distribution function as 𝜙𝑦=1√2𝜋𝜈²exp−𝑦−𝑚²2𝜈²,𝑦∈ℝ, (2) which provides an important averaging tool for the unobserved process 𝑌_𝑡 as documented well in [13]. Notation ⟨⋅⟩ is going to be used for the expectation with respect to this invariant distribution.

Based on Lando [5], the price of zero-recovery defaultable bond is given by the reduced form 𝑃^𝜖,𝛿𝑡,𝑟,𝜆,𝑦,𝑧;𝑇=𝐸^∗𝑒^{−∫𝑇𝑡(𝑟𝑠+𝜆𝑠)𝑑𝑠}∣𝑟_𝑡=𝑟,𝜆_𝑡=𝜆,𝑌_𝑡=𝑦,𝑍_𝑡=𝑧 (3) under the risk-neutral probability measure. Here, the Markov property of the joint process (𝑟_𝑡,𝜆_𝑡,𝑌_𝑡,𝑍_𝑡) was used. Of course, the zero-recovery assumption is not appropriate from an economic point of view but, for the purpose of mathematical simplicity, we assume in this paper that the loss rate equals one identically. Then, using the four-dimensional Feynman-Kac formula (cf. [16]), we obtain a singularly and regularly perturbed partial differential equation (PDE) problem given by ℒ^𝜖,𝛿𝑃^𝜖,𝛿ℒ𝑡,𝑟,𝜆,𝑦,𝑧;𝑇=0,𝑡<𝑇,_𝜖,𝛿1∶=𝜖ℒ₀+1√𝜖ℒ₁+ℒ₂+√𝛿𝒦₁+𝛿𝒦₂+𝛿𝜖𝒦₃,𝑃^𝜖,𝛿𝑡,𝑟,𝜆,𝑦,𝑧;𝑇_𝑡=𝑇=1, (4) where ℒ₀=𝜕𝑚−𝑦𝜕𝑦+𝜈²𝜕²𝜕𝑦²,ℒ₁=√2𝜎𝜈𝜌_𝑟𝑦√𝑟𝜕²+√𝜕𝑟𝜕𝑦2𝑓𝑦,𝑧𝜈𝜌_𝜆𝑦√𝜆𝜕²,ℒ𝜕𝜆𝜕𝑦₂=𝜕𝑟𝜕𝑡+𝑎^∗𝜕−𝑟+1𝜕𝑟2𝜎²𝑟𝜕²𝜕𝑟²+𝑎𝜆^∗𝜕−𝜆+1𝜕𝜆2𝑓²𝜆𝜕𝑦,𝑧²𝜕𝜆²−()𝒦𝑟+𝜆𝐼,₁=𝜎ℎ(𝑧)𝜌_𝑟𝑧√𝑟𝜕²ℎ(𝑧)𝜌𝜕𝑟𝜕𝑧+𝑓𝑦,𝑧_𝜆𝑧√𝜆𝜕²,𝒦𝜕𝜆𝜕𝑧₂(𝑧)𝜕=𝑔+1𝜕𝑧2ℎ²(𝑧)𝜕²𝜕𝑧²,𝒦₃=√2𝜈ℎ(𝑧)𝜌_𝑦𝑧𝜕².𝜕𝑦𝜕𝑧 (5) Note that (1/𝜖)ℒ₀ is the infinitesimal generator of the OU process 𝑌_𝑡. The operator ℒ₁ contains the mixed partial derivative due to the correlation between 𝑟_𝑡 and 𝑌_𝑡 as well as between 𝜆_𝑡 and 𝑌_𝑡. ℒ₂ is the operator of the canonical two-factor CIR model with volatility at the volatility level 𝑓(𝑦,𝑧). 𝒦₁ includes the mixed partial derivatives due to the correlation between 𝑟_𝑡 and 𝑍_𝑡 and between 𝜆_𝑡 and 𝑍_𝑡. 𝒦₂ is the infinitesimal generator of the process 𝑍_𝑡. 𝒦₃ holds the mixed partial derivative due to the correlation between 𝑌_𝑡 and 𝑍_𝑡.

3. Multiscale Analysis

Since it is difficult to solve the PDE problem (4) itself, we are interested in the following asymptotic expansions: 𝑃^𝜖,𝛿𝑡,𝑟,𝜆,𝑦,𝑧;𝑇=^∞_𝑗=0𝛿^𝑗/2𝑃^𝜖_𝑗,𝑃𝑡,𝑟,𝜆,𝑦,𝑧;𝑇_𝜖𝑗𝑡,𝑟,𝜆,𝑦,𝑧;𝑇=^∞_𝑖=0𝜖^𝑖/2𝑃_𝑖,𝑗𝑡,𝑟,𝜆,𝑦,𝑧;𝑇, (6)(7) so that 𝑃^𝜖,𝛿 is a series of the general term 𝜖^𝑖/2𝛿^𝑗/2𝑃_𝑖,𝑗. Plugging the expansion (6) into (4) leads to 𝑃^𝜖₀ and 𝑃^𝜖₁ given by the solutions of the PDEs 1𝜖ℒ₀+1√𝜖ℒ₁+ℒ₂𝑃^𝜖₀𝑃𝑡,𝑟,𝜆,𝑦,𝑧;𝑇=0,𝑡<𝑇,_𝜖0(𝑡,𝑟,𝜆,𝑦,𝑧;𝑇)_𝑡=𝑇1=1,𝜖ℒ₀+1√𝜖ℒ₁+ℒ₂𝑃^𝜖₁𝒦𝑡,𝑟,𝜆,𝑦,𝑧;𝑇=−₁+1√𝜖𝒦₃𝑃^𝜖₀𝑃,𝑡<𝑇,_𝜖1(𝑡,𝑟,𝜆,𝑦,𝑧;𝑇)_𝑡=𝑇=0, (8)(9) respectively. From now on, we employ an analytic technique of [13] to approximate the two terms 𝑃^𝜖₀ and 𝑃^𝜖₁ based on the ergodic property of the Ornstein-Uhlenbeck (OU) process 𝑌_𝑡.

3.1. The Leading Order Term

Applying the expansion (7) with 𝑗=0 to (8) leads to 1𝜖ℒ₀𝑃_0,0+1√𝜖ℒ₀𝑃_1,0+ℒ₁𝑃_0,0+ℒ₀𝑃_2,0+ℒ₁𝑃_1,0+ℒ₂𝑃_0,0+√𝜖ℒ₀𝑃_3,0+ℒ₁𝑃_2,0+ℒ₂𝑃_1,0+⋯=0 (10) from which we obtain an affine form of the leading order term 𝑃_0,0 as follows.

Theorem 1. Assume that the partial derivative of 𝑃_𝑖,𝑗 with respect to 𝑦 does not grow as much as 𝜕𝑃_𝑖,𝑗/𝜕𝑦∼𝑒^𝑦2/2 as 𝑦 goes to infinity. Then the leading order term 𝑃₀∶=𝑃_0,0 of the expansion (7) with 𝑗=0 is independent of the fast scale variable 𝑦 and further it has the affine representation 𝑃₀(𝑇−𝜏,𝑟,𝜆,𝑧;𝑇)=𝐶̆(𝜏,𝑧)𝑒^{̆̆−𝑟𝐴(𝜏)−𝜆𝐵(𝜏)} (11) with 𝜏=𝑇−𝑡, 𝐴̆(0)=𝐵̆(0)=0 and 𝐶̆(0,𝑧)=1, where 𝐴̆(𝜏), 𝐵̆(𝜏), and 𝐶̆(𝜏,𝑧) are given by, 𝐴̆(𝜏)=2𝑒^𝛾1𝜏−12𝛾₁+𝛾₁+𝑎(𝑒^𝛾1𝜏),𝐵̆(𝜏)=2𝑒−1_𝛾2𝜏−12𝛾₂+𝛾₂+𝑎(𝑒^𝛾2𝜏),𝐶̆()=−1𝜏,𝑧2𝛾₁𝑒^{(1/2)(𝛾1+𝑎)𝜏}2𝛾₁+𝛾₁+𝑎(𝑒^𝛾1𝜏−1)^{2𝑎𝑟∗/𝜎2}×2𝛾₂𝑒^(1/2)(𝛾2+𝑎_)𝜏2𝛾₂+(𝛾₂+𝑎)(𝑒^𝛾2𝜏−1)²𝑎_𝜆^∗_/^̆_𝜎²_(𝑧),𝛾₁∶=√𝑎²+2𝜎²,𝛾₂∶=𝑎²+2𝜎̆²(𝑧), (12) respectively. Here, 𝜎̆ is defined by 𝜎̆²(𝑧)∶=⟨𝑓²⟩(𝑧)=_ℝ𝑓²𝑦,𝑧𝜙𝑦𝑑𝑦 (13) in terms of 𝜙 (the invariant distribution of 𝑌_𝑡).

Proof. Multiplying (10) by 𝜖 and letting 𝜖 go to zero, we obtain the ordinary differential equation (ODE) ℒ₀𝑃₀=0. Recalling that the operator ℒ₀ is the generator of the OU process 𝑌_𝑡, the solution 𝑃₀ of this ODE must be independent of the 𝑦 variable due to the assumed growth condition; 𝑃₀=𝑃₀(𝑡,𝑟,𝜆,𝑧;𝑇). From the 𝑂(1/√𝜖) terms of (10), we have the PDE ℒ₀𝑃_1,0+ℒ₁𝑃₀=0. If we apply the 𝑦-independence of 𝑃₀ to this PDE, then it reduces to the ODE ℒ₀𝑃_1,0=0 and so 𝑃_1,0 is independent of 𝑦 by the same reason as in the case of 𝑃₀; 𝑃_1,0=𝑃_1,0(𝑡,𝑟,𝜆,𝑧;𝑇). So far, the two terms 𝑃₀ and 𝑃_1,0 do not depend on the current level 𝑦 of the fast scale volatility driving the process 𝑌_𝑡.
One can continue to eliminate the terms of order 1, √𝜖, 𝜖, and so forth. From the 𝑂(1) terms of (10), we get ℒ₀𝑃_2,0+ℒ₁𝑃_1,0+ℒ₂𝑃₀=0. This PDE becomes ℒ₀𝑃_2,0+ℒ₂𝑃₀=0 (14) due to the 𝑦-independence of 𝑃_1,0 obtained above. Equation (14) is a Poisson equation for 𝑃_2,0 with respect to the operator ℒ₀ in the 𝑦 variable. It is well known (cf. [17]) from the Fredholm alternative (solvability condition) that it has a solution only if ℒ₂𝑃₀ is centered with respect to the invariant distribution of 𝑌_𝑡; that is, ⟨ℒ₂⟩𝑃₀ℒ=0,₂𝜕∶=𝑟𝜕𝑡+𝑎^∗𝜕−𝑟+1𝜕𝑟2𝜎²𝑟𝜕²𝜕𝑟²+𝑎𝜆^∗𝜕−𝜆+1𝜕𝜆2𝜎̆²(𝑧)𝜆𝜕²𝜕𝜆²−(𝑟+𝜆)𝐼. (15) Noting that (15) is nothing but a PDE corresponding to the intensity with constant volatility replaced by the function 𝜎̆(𝑧), the affine solution of (15) is given by (11)-(12) from the well-known solution in [6].

3.2. The Correction Terms

Next, we derive the first-order correction terms 𝑃_1,0 and 𝑃_0,1 from the leading order solution 𝑃₀ obtained above. Contrary to the result obtained in [12], the correction terms are not given by closed form solutions in the present case of a CIR model (instead of a Vasicek model) for the interest rate.

Theorem 2. Assume that the partial derivative of 𝑃_𝑖,𝑗 with respect to 𝑦 does not grow as much as 𝜕𝑃_𝑖,𝑗/𝜕𝑦∼𝑒^𝑦2/2 as 𝑦 goes to infinity. The correction term 𝑃_1,0 is independent of the variable 𝑦 and 𝑃^𝜖_1,0∶=√𝜖𝑃_1,0(𝑡,𝑟,𝜆,𝑧;𝑇) is given by 𝑃^𝜖_1,0(𝑇−𝜏,𝑟,𝜆,𝑧;𝑇)=𝐷̆(𝜏,𝑟,𝜆,𝑧)𝑃₀(𝑇−𝜏,𝑟,𝜆,𝑧;𝑇) (16) with 𝐷̆(0,𝑟,𝜆,𝑧)=0, where 𝐷̆(𝜏,𝑟,𝜆,𝑧) is given by the solution of the PDE 𝜕𝐷̆−𝑎𝑟𝜕𝜏^∗−𝑟−𝜎²𝑟𝐴̆(𝜏𝜕𝐷̆)−1𝜕𝑟2𝜎²𝑟𝜕²𝐷̆𝜕𝑟²−𝑎𝜆^∗−𝜆−𝜎̆²(𝑧)𝜆𝐵̆(𝜏𝜕𝐷̆)−1𝜕𝜆2𝜎̆²(𝑧)𝜆𝜕²𝐷̆𝜕𝜆²(),𝑙()(𝑧)√=𝑙𝜏,𝑟,𝜆,𝑧𝜏,𝑟,𝜆,𝑧∶=𝑞𝜆𝐵̆³(𝜏)−𝑠(𝑧)√𝑟𝐴̆(𝜏)𝐵̆²(𝜏),𝑞(𝑧)𝜈√∶=𝜖√2𝜌_𝜆𝑦√𝜆⟨𝑓𝜃_𝑦𝑠(𝑧)𝜈√⟩,∶=𝜖√2𝜌_𝑟𝑦𝜎√𝑟⟨𝜃_𝑦⟩. (17) Here, 𝐴̆(𝜏) and 𝐵̆(𝜏) are given by (12) and the function 𝜃∶ℝ²→ℝ is defined by the solution of ℒ₀𝜃𝑦,𝑧=𝑓²𝑦,𝑧−⟨𝑓²⟩(𝑧) (18) and 𝜃_𝑦 denotes the partial derivative with respect to the 𝑦 variable.

Proof. The 𝑂(√𝜖) terms of (10) lead to ℒ₀𝑃_3,0+ℒ₁𝑃_2,0+ℒ₂𝑃_1,0=0 which is a Poisson equation for 𝑃_3,0 and so the solvability condition leads to the PDE ℒ₁𝑃_2,0+ℒ₂𝑃_1,0=0. (19) Meanwhile, from (14) and (15), we have 𝑃_2,0=−ℒ₀⁻¹ℒ₂−⟨ℒ₂⟩𝑃₀+𝑐(𝑡,𝑟,𝜆,𝑧) (20) for some function 𝑐(𝑡,𝑟,𝜆,𝑧) that does not rely on the 𝑦 variable. Plugging (20) into (19), a PDE for 𝑃_1,0 is obtained as ⟨ℒ₂⟩𝑃_1,0=⟨ℒ₁ℒ₀⁻¹ℒ₂−⟨ℒ₂⟩⟩𝑃₀ (21) with the terminal condition 𝑃_1,0(𝑡,𝑟,𝜆,𝑧;𝑇)|_𝑡=𝑇=0. This implies that 𝑃_1,0 is 𝑦-independent. Since we focus on the first-order correction terms to 𝑃₀, we reset (21) with respect to 𝑃^𝜖_1,0 leading to ⟨ℒ₂⟩𝑃^𝜖_1,0=ℋ₁𝑃₀,ℋ₁∶=√𝜖ℒ₁ℒ₀⁻¹ℒ₂−ℒ₂. (22) To calculate the operator ℋ₁, we first find ℒ₂−⟨ℒ₂1⟩=2𝑓²𝑦,𝑧−⟨𝑓²⟩(𝑧)𝜕²𝜕𝜆² (23) from (5) and so the operator ℋ₁ is expressed as ℋ₁=𝑞(𝑧)√𝜆𝜕³_𝜆𝜆𝜆+𝑠(𝑧)√𝑟𝜕³_𝑟𝜆𝜆, where 𝑞(𝑧) and 𝑠(𝑧) are the ones given in the theorem. Then, from (15), Theorem 1 and the change of variables 𝜏=𝑇−𝑡, the PDE (22) leads to 𝜕𝑃^𝜖_1,0𝑟𝜕𝜏=𝑎^∗𝜕𝑃−𝑟^𝜖1,0+1𝜕𝑟2𝜎²𝑟𝜕²𝑃^𝜖_1,0𝜕𝑟²+𝑎𝜆^∗𝜕𝑃−𝜆^𝜖1,0+1𝜕𝜆2𝜎̆²(𝑧)𝜆𝜕²𝑃^𝜖_1,0𝜕𝜆²−(𝑟+𝜆)𝑃^𝜖_1,0+𝐶̆(𝜏,𝑧)𝑒^{̆̆−𝑟𝐴(𝜏)−𝜆𝐵(𝜏)}×𝑞(𝑧)√𝜆𝐵̆³(𝜏)+𝑠(𝑧)√𝑟𝐴̆(𝜏)𝐵̆²(𝜏) (24) with the initial condition 𝑃^𝜖_1,0(𝑇−𝜏,𝑟,𝜆,𝑧;𝑇)|_𝜏=0=0. Finally, plugging 𝑃^𝜖_1,0 of the form (16) into (24), we obtain the result of Theorem 2 by direct computation.

Similarly, one can derive the correction term 𝑃_0,1 also as follows.

Theorem 3. Assume that the partial derivative of 𝑃_𝑖,𝑗 with respect to 𝑦 does not grow as much as 𝜕𝑃_𝑖,𝑗/𝜕𝑦∼𝑒^𝑦2/2 as 𝑦 goes to infinity. The first correction term 𝑃_0,1 does not depend on the variable 𝑦 and 𝑃^𝛿_0,1∶=√𝛿𝑃_0,1(𝑡,𝑟,𝜆,𝑧;𝑇) is given by 𝑃^𝛿_0,1(𝑇−𝜏,𝑟,𝜆,𝑧;𝑇)=𝐸̆(𝜏,𝑟,𝜆,𝑧)𝑃₀(𝑇−𝜏,𝑟,𝜆,𝑧;𝑇) (25) with 𝐸̆(0,𝑟,𝜆,𝑧)=0, where 𝐸̆(𝜏,𝑟,𝜆,𝑧) is given by the solution of the PDE 𝜕𝐸̆−𝑎𝑟𝜕𝜏^∗−𝑟−𝜎²𝑟𝐴̆(𝜏𝜕𝐸̆)−1𝜕𝑟2𝜎²𝑟𝜕²𝐸̆𝜕𝑟²−𝑎𝜆^∗−𝜆−𝜎̆²(𝑧)𝜆𝐵̆(𝜏𝜕𝐸̆)−1𝜕𝜆2𝜎̆²(𝑧)𝜆𝜕²𝐸̆𝜕𝜆²(),𝑚()1=𝑚𝜏,𝑟,𝜆,𝑧𝜏,𝑟,𝜆,𝑧∶=𝐶̆𝑢(𝑧)√𝑟𝐴̆(𝜏)𝜕𝐶̆(𝑧)√𝜕𝑧+𝑣𝜆𝐵̆(𝜏)𝜕𝐶̆,𝑢(𝑧)√𝜕𝑧∶=−𝛿𝜌_𝑟𝑧(𝑧),𝑣(𝑧)√𝜎ℎ∶=−𝛿𝜌_𝜆𝑧𝜎(𝑧)ℎ(𝑧). (26) Here, 𝐴̆(𝜏), 𝐵̆(𝜏), and 𝐶̆(𝜏,𝑧) are given by (12) and 𝜎 is defined by 𝜎(𝑧)∶=⟨𝑓⟩(𝑧)=_ℝ𝑓𝑦,𝑧𝜙𝑦𝑑𝑦 (27) in terms of 𝜙 (the invariant distribution of 𝑌_𝑡).

Proof. Applying the expansion (7) with 𝑗=0 and 𝑗=1 to (9) results in 1𝜖ℒ₀𝑃_0,1+1√𝜖ℒ₀𝑃_1,1+ℒ₁𝑃_0,1+ℒ₀𝑃_2,1+ℒ₁𝑃_1,1+ℒ₂𝑃_0,1+√𝜖ℒ₀𝑃_3,1+ℒ₁𝑃_2,1+ℒ₂𝑃_1,11+⋯=−√𝜖𝒦₃𝑃₀−𝒦₁𝑃₀+𝒦₃𝑃_1,0−√𝜖𝒦₁𝑃_1,0+𝒦₃𝑃_2,0+⋯. (28) Then, from the 𝑂(1/𝜖) term of (28), 𝑃_0,1 is 𝑦-independent as the solution of the ODE ℒ₀𝑃_0,1=0 under the assumed growth condition. Since both 𝑃₀ (obtained in Theorem 1) and 𝑃_0,1 are 𝑦-independent, from the 𝑂(1/√𝜖) terms of (28), we have the ODE ℒ₀𝑃_1,1=0. So, by the same reason as in the case of 𝑃_0,1=0, 𝑃_1,1 is also independent of 𝑦. Hence, the two terms 𝑃_0,1 and 𝑃_1,1 do not depend on the current level 𝑦 of the fast scale volatility of intensity.
One can continue to remove the terms of order 1, √𝜖, 𝜖, and so forth. From the 𝑂(1) terms and the 𝑦-independence of 𝑃_1,0 and 𝑃_1,1, we have the PDE ℒ₀𝑃_2,1+ℒ₂𝑃_0,1+𝒦₁𝑃₀=0. This is a Poisson equation for 𝑃_2,1 with respect to the operator ℒ₀ in the 𝑦 variable and so it has a solution only if ⟨ℒ₂⟩𝑃_0,1=−⟨𝒦₁⟩𝑃₀ holds. If we reset this PDE with respect to 𝑃^𝛿_0,1, then we have ℒ₂𝑃^𝛿_0,1=ℋ₂𝑃₀,ℋ₂∶=−√𝛿𝒦₁, (29) where the operator ℋ₂ is the same as ℋ₂=𝑢(𝑧)√𝑟𝜕²_𝑟𝑧+𝑣(𝑧)√𝜆𝜕²_𝜆𝑧 from the definition of 𝒦₁ given by (5).
Using the operator (15), Theorem 1, and the change of variable 𝜏=𝑇−𝑡, we write the PDE (29) as follows: 𝜕𝑃^𝛿_0,1𝑟𝜕𝜏=𝑎^∗𝜕𝑃−𝑟^𝛿0,1+1𝜕𝑟2𝜎²𝑟𝜕²𝑃^𝛿_0,1𝜕𝑟²+𝑎𝜆^∗𝜕𝑃−𝜆^𝛿0,1+1𝜕𝜆2𝜎̆²(𝑧)𝜆𝜕²𝑃^𝛿_0,1𝜕𝜆²−(𝑟+𝜆)𝑃^𝛿_0,1+𝑒^{̆̆−𝑟𝐴(𝜏)−𝜆𝐵(𝜏)}𝑢(𝑧)√𝑟𝐴̆(𝜏)𝜕𝐶̆(𝑧)√𝜕𝑧+𝑣𝜆𝐵̆(𝜏)𝜕𝐶̆,𝜕𝑧 (30) where the initial condition 𝑃^𝛿_0,1(𝑇−𝜏,𝑟,𝜆,𝑧;𝑇)|_𝜏=0=0 is satisfied. Plugging the form (25) into (30), we obtain the result of Theorem 3 by direct computation.

3.3. Accuracy

Synthesizing Theorems 1, 2, and 3, we obtain the following asymptotic representation of the price of the defaultable bond at time 𝑡: 𝑃^𝜖,𝛿𝑡,𝑟,𝜆,𝑦,𝑧;𝑇≈𝑃^𝜖,𝛿∶=𝑃₀+𝑃^𝜖_1,0+𝑃^𝛿_0,1=𝐷̆()+𝐸̆()×𝐶̆()𝑒1+𝜏,𝑟,𝜆,𝑧𝜏,𝑟,𝜆,𝑧𝜏,𝑧^{̆̆−𝑟𝐴(𝜏)−𝜆𝐵(𝜏)}, (31) where 𝐴̆(𝜏), 𝐵̆(𝜏), 𝐶̆(𝜏,𝑧), 𝐷̆(𝜏,𝑟,𝜆,𝑧), and 𝐸̆(𝜏,𝑟,𝜆,𝑧) are given by (12), (17), and (26), respectively. We note that the function 𝑔(𝑧) (the drift term of the slow scale variation) in (1) does not affect the present form of the approximated bond price 𝑃^𝜖,𝛿 due to the order of 𝛿 in front of 𝑔(𝑧) in (1). It should influence the next order approximation.

In this approximation, the leading order price is determined by the functions 𝐴̆(𝜏), 𝐵̆(𝜏), and 𝐶̆(𝜏,𝑧) and the first-order corrections are given by 𝐷̆(𝜏,𝑟,𝜆,𝑧) and 𝐸̆(𝜏,𝑟,𝜆,𝑧). Here, 𝐷̆(𝜏,𝑟,𝜆,𝑧) is determined by the two group parameters 𝑞(𝑧) and 𝑠(𝑧) and 𝐸̆(𝜏,𝑟,𝜆,𝑧) is determined by 𝑢(𝑧) and 𝑣(𝑧). These four group parameters absorb the original parameters and functions related to 𝑌_𝑡 and 𝑍_𝑡. The group parameter 𝑞(𝑧) absorbs 𝜖, 𝜈, 𝜌_𝜆𝑦, and 𝑓 while 𝑠(𝑧) absorbs 𝜖, 𝜈, 𝜌_𝑟𝑦, and 𝑓. The group parameter 𝑢(𝑧) absorbs 𝛿, 𝜌_𝑟𝑧, and ℎ while 𝑣(𝑧) absorbs 𝛿, 𝜌_𝜆𝑧, and ℎ. Particularly, the functions 𝑓, 𝑔 and ℎ determining the original model (1) need not be specified to price the defaultable bond. The asymptotic expansions based upon the ergodic property of the OU process 𝑌_𝑡 provide the reduction of the model parameters and the absorption of the model functions into the group parameters.

All the derivations given so far are formal (as usually done in this kind of research work). So, we discuss the accuracy of the asymptotic approximation (31) as shown in the following theorem.

Theorem 4. Let 𝑃^𝜖,𝛿 be defined by (31). Then for any fixed (𝑡,𝑟,𝜆,𝑦,𝑧) there exists a positive constant 𝐶, which depends on (𝑡,𝑟,𝜆,𝑦,𝑧) but not on 𝜖 and 𝛿, such that 𝑃^𝜖,𝛿𝑡,𝑟,𝜆,𝑦,𝑧;𝑇−𝑃^𝜖,𝛿(𝑡,𝑟,𝜆,𝑧;𝑇)≤𝐶(𝜖+𝛿) (32) holds for all 0<𝜖,  𝛿≤1.

Proof. To prove the inequality (32), we first define a residual 𝑅^𝜖,𝛿 by writing 𝑃^𝜖,𝛿 as 𝑃^𝜖,𝛿=^𝑖,𝑗=3_𝑖,𝑗=0𝜖^𝑖/2𝛿^𝑗/2𝑃_𝑖,𝑗−𝑅^𝜖,𝛿. (33) The second step is to compute ℒ^𝜖,𝛿𝑅^𝜖,𝛿 using the properties of 𝑃^𝜖,𝛿,  𝑃₀,…,𝑃_3,3 obtained above. The last step is to write 𝑅^𝜖,𝛿 as a probabilistic representation (the Feynman-Kac formula) and to show the desired estimate (32). This type of argument is standard and our formulation would not make a different impact on the derivation in Theorem 3.1 of [12]. So, we leave the details there.

4. Numerical Experiment

In this section, we compute numerically the fast term factor 𝐷̆ and the slow term factor 𝐸̆ of the stochastic volatility of default intensity and investigate the impact of the multiscale factors on the functional behavior of the bond price with the constant volatility of default intensity.

Table 1 shows the average cumulative default rates of bond issuers for the period of 1983–2009. It says, for instance, that over a ten-year 𝐵-rated issuers were in default at a 44.982% average rate between 1983 and 2009. Using this table, in general, one can calculate the average default intensity between time 0 and 𝑡 from the cumulative probability of default by time 𝑡. For example, the average 20-year default intensity for 𝐵-rated issuers is given by −(1/20)ln(1−0.65493)=0.0532.

Table 1 

Average cumulative default rates (%), 1983–2009 (source: Moody's).

Based on Table 1 and the data [13], we choose the parameter values 𝑎=1.2,  𝑎=1.3,  𝜎=0.06,  𝜎̆²=0.03,  𝑟^∗=0.1,  𝜆=0.0455,  and  𝜆^∗=0.0532. Maturity runs from 0 to 20 years and the interest rate is fixed as 𝑟=0.07. Figures 1 and 2 show the bond price correction term structure depending on the group parameters (𝑞,𝑠) and (𝑢,𝑣), respectively, for the fast term factor 𝐷̆ and the slow term factor 𝐸̆. One can notice that the fast scale correction creates a more significant impact on the bond price than the slow scale correction and the impact tends to increase as time to maturity increases.

Figure 1 

The price correction term structure factor 𝐷̆ is shown by dotted, dashed, and sold lines as functions of time to maturity, respectively, for (𝑞,𝑠)=(0.2,0.1), (0.3,0.4), and (0.6,0.7).

Figure 2 

The price correction term structure factor 𝐸̆ is shown by dotted, dashed and sold lines as functions of time to maturity, respectively, for (𝑢,𝑣)=(0.1,0.2), (0.3,0.5), and (0.7,0.8).

5. Bond Option Pricing

In this section, we are interested in an asymptotic pricing formula for a European option, where the underlying asset itself is a zero-coupon bond with default risk.

We use notations 𝑇 and 𝑇₀,  𝑇₀<𝑇, to denote the maturity of the defaultable bond and the maturity of a European option written on the defaultable bond, respectively. It is assumed that the option becomes invalid when a default occurs prior to 𝑇₀. The short rate process 𝑟_𝑡, the intensity process 𝜆_𝑡, the fast volatility process 𝑌_𝑡 of the intensity, and the slow volatility process 𝑍_𝑡 of the intensity are given by the SDEs (1).

The option price at time 𝑡 for an observed short rate 𝑟_𝑡=𝑟, an intensity level 𝜆_𝑡=𝜆, a fast volatility level 𝑌_𝑡=𝑦, and a slow volatility level 𝑍_𝑡=𝑧, denoted by 𝑄(𝑡,𝑟,𝜆,𝑦,𝑧;𝑇₀,𝑇), is defined by 𝑄𝑡,𝑟,𝜆,𝑦,𝑧;𝑇₀,𝑇=𝐸^∗𝑒^{−∫𝑇0𝑡𝑟𝑠𝑑𝑠}ℎℬ𝑇₀,𝑇∣𝑟_𝑡=𝑟,𝜆_𝑡=𝜆,𝑌_𝑡=𝑦,𝑍_𝑡=𝑧 (34) under an equivalent martingale measure, where the bond price ℬ(𝑇₀,𝑇) is given by ℬ𝑇₀,𝑇=𝑃^𝜖,𝛿𝑇₀,𝑟_𝑇0,𝜆_𝑇0,𝑌_𝑇0,𝑍_𝑇0;𝑇=𝐸^∗𝑒^{−∫𝑇𝑇0(𝑟𝑠+𝜆𝑠)𝑑𝑠}∣𝑟_𝑇0,𝜆_𝑇0,𝑌_𝑇0,𝑍_𝑇0 (35) and ℎ is the payoff function of the option at time 𝑇₀. Here, it is assumed that the payoff function ℎ is at best linearly growing at infinity and is a smooth function. This smoothness assumption may be too severe in practical point of view as ℎ is not differentiable at the exercise price in the case of classical European call or put option. The smoothness assumption on ℎ, however, can be removed similarly to the argument in [18] or [12]. So, for simplicity, ℎ is assumed to be a smooth function.

By the four-dimensional Feynman-Kac formula, we transform the above integral form (34) into the PDE 𝜕𝑄𝑟𝜕𝑡+𝑎^∗−𝑟𝜕𝑄+1𝜕𝑟2𝜎²𝑟𝜕²𝑄𝜕𝑟²+𝑎𝜆^∗−𝜆𝜕𝑄+1𝜕𝜆2𝑓²𝜆𝜕𝑦,𝑧²𝑄𝜕𝜆²+1𝜖𝑚−𝑦𝜕𝑄+𝜈𝜕𝑦𝜖²𝜕²𝑄𝜕𝑦²+𝛿𝑔(𝑧)𝜕𝑄+1𝜕𝑧2𝛿ℎ²(𝑧)𝜕²𝑄𝜕𝑧²+𝜌_𝑟𝑦𝜎√𝑟𝜈√2√𝜖𝜕²𝑄𝜕𝑟𝜕𝑦+𝜌_𝑟𝑧𝜎√𝑟√(𝑧)𝜕𝛿ℎ²𝑄𝜕𝑟𝜕𝑧+𝜌_𝜆𝑦𝑓𝑦,𝑧√𝜆𝜈√2√𝜖𝜕²𝑄𝜕𝜆𝜕𝑦+𝜌_𝜆𝑧𝑓𝑦,𝑧√𝜆√(𝑧)𝜕𝛿ℎ²𝑄𝜕𝜆𝜕𝑧+𝜌_𝑦𝑧𝜈√2√𝜖√𝛿ℎ(𝑧)𝜕²𝑄−()𝜕𝑦𝜕𝑧𝑟+𝜆𝑄=0 (36) with the terminal condition 𝑄(𝑡,𝑟,𝜆,𝑦,𝑧;𝑇₀,𝑇)|_𝑡=𝑇0=ℎ(𝑃^𝜖,𝛿(𝑇₀,𝑟,𝜆,𝑦,𝑧;𝑇)). Then, keeping the notation used in Section 3 but with a different terminal condition, the option price 𝑄(𝑡,𝑟,𝜆,𝑦,𝑧;𝑇₀,𝑇) is given by the solution of the PDE problem ℒ^𝜖,𝛿𝑄𝑡,𝑟,𝜆,𝑦,𝑧;𝑇₀𝑄,𝑇=0,𝑡,𝑟,𝜆,𝑦,𝑧;𝑇₀,𝑇_𝑡=𝑇0=ℎ𝑃^𝜖,𝛿𝑇₀,𝑟,𝜆,𝑦,𝑧;𝑇. (37) Taking the asymptotic expansions 𝑄=𝑄^𝜖₀+√𝛿𝑄^𝜖₁+𝛿𝑄^𝜖₂𝑄+⋯,_𝜖𝑘=𝑄_0,𝑘+√𝜖𝑄_1,𝑘+𝜖𝑄_2,𝑘+𝜖^3/2𝑄_3,𝑘𝑄+⋯,𝑘=0,1,2,…,𝑡,𝑟,𝜆,𝑦,𝑧;𝑇₀,𝑇_𝑡=𝑇0=ℎ𝑃₀𝑇₀+𝑃,𝑟,𝜆,𝑧;𝑇_𝜖1,0𝑇₀,𝑟,𝜆,𝑧;𝑇ℎ^′𝑃₀𝑇₀+𝑃,𝑟,𝜆,𝑧;𝑇_𝛿0,1𝑇₀,𝑟,𝜆,𝑧;𝑇ℎ^′𝑃₀𝑇₀,𝑟,𝜆,𝑧;𝑇+⋯, (38)(39)(40) where 𝑃₀, 𝑃^𝜖_1,0, and 𝑃^𝛿_0,1 are given by Theorems 1, 2, and 3, respectively, we employ the same asymptotic analysis as in Section 3 to derive an asymptotic pricing formula for the bond option.

First, the terms of order 1/𝜖 and 1/√𝜖 in the asymptotic PDE (37) yield the 𝑦-independence of 𝑄_0,0 (in brief, 𝑄₀), 𝑄_1,0, and 𝑄_0,1 under the same growth condition as in Section 3.

The 𝑂(1) terms in (37) give a Poisson equation for 𝑄_2,0 from which the solvability condition ⟨ℒ₂⟩𝑄₀=0 has to be satisfied. From (40), the corresponding terminal condition is given by 𝑄₀|_𝑡=𝑇0=ℎ(𝑃₀(𝑇₀,𝑟,𝜆,𝑧;𝑇)). Then we have 𝑄₀𝑡,𝑟,𝜆,𝑧;𝑇₀,𝑇=𝐸^∗𝑒^{−∫𝑇0𝑡𝑟𝑠𝑑𝑠}ℎ𝑃₀𝑇₀,𝑟,𝜆,𝑧;𝑇∣𝑟_𝑡=𝑟,𝜆_𝑡=𝜆,𝑍_𝑡=𝑧, (41) where 𝑃₀(𝑇₀,𝑟,𝜆,𝑧;𝑇) is given by (11) in Theorem 1 at time 𝑡=𝑇₀.

The 𝑂(√𝜖) terms in (37) lead to a Poisson equation for 𝑄_3,0, where the solvability condition is given by ⟨ℒ₁𝑄_2,0+ℒ₂𝑄_1,0⟩=0. If we put 𝑄_1,0∶=√𝜖𝑄_1,0, then this solvability condition leads to ℒ₂𝑄_1,0𝑡,𝑟,𝜆,𝑧;𝑇₀,𝑇=ℋ₁𝑄₀𝑡,𝑟,𝜆,𝑧;𝑇₀,𝑇 (42) with a terminal condition given by 𝑄_1,0𝑡,𝑟,𝜆,𝑧;𝑇₀,𝑇|_𝑡=𝑇0=𝑃^𝜖_1,0𝑇₀,𝑟,𝜆,𝑧;𝑇ℎ^′𝑃₀𝑇₀,𝑟,𝜆,𝑧;𝑇 (43) from (40), where the operator ℋ₁ is given by (22) and 𝑃^𝜖_1,0 is given by the solution (16) in Theorem 2. Then, by applying the Feynman-Kac formula to (42) and (43), we obtain the following probabilistic representation: 𝑄_1,0𝑡,𝑟,𝜆,𝑧;𝑇₀,𝑇=𝐸^∗𝑒^{−∫𝑇0𝑡𝑟𝑠𝑑𝑠}𝑃^𝜖_1,0𝑇₀,𝑟,𝜆,𝑧;𝑇ℎ^′𝑃₀𝑇₀,𝑟,𝜆,𝑧;𝑇−^𝑇0_𝑡𝑒^{−∫𝑢𝑡𝑟𝑠𝑑𝑠}ℋ₁𝑄₀𝑢,𝑟,𝜆,𝑧;𝑇₀,𝑇𝑑𝑢∣𝑟_𝑡𝜆=𝑟,_𝑡=𝜆,𝑍_𝑡=𝑧. (44) Now, the 𝑂(√𝛿) terms in (37) lead to a Poisson equation for 𝑄_2,1 such that the solvability condition ⟨ℒ₂𝑄_0,1⟩=−⟨𝒦₁𝑄₀⟩ holds. If we let 𝑄_0,1∶=√𝛿𝑄_0,1, then this solvability condition leads to ⟨ℒ₂⟩𝑄_0,1𝑡,𝑟,𝜆,𝑧;𝑇₀,𝑇=ℋ₂𝑄₀𝑡,𝑟,𝜆,𝑧;𝑇₀,𝑇, (45) where the operator ℋ₂ is given by (29). Also, from the terminal condition (40), we have 𝑄_0,1𝑡,𝑟,𝜆,𝑧;𝑇₀,𝑇_𝑡=𝑇0=𝑃^𝛿_0,1𝑇₀,𝑟,𝜆,𝑧;𝑇ℎ^′𝑃₀𝑇₀,𝑟,𝜆,𝑧;𝑇, (46) where 𝑃^𝛿_0,1 is given by the solution (25) in Theorem 3. By applying the Feynman-Kac formula to (45) and (46), we have the following probabilistic representation: 𝑄_0,1𝑡,𝑟,𝜆,𝑧;𝑇₀,𝑇=𝐸^∗𝑒^{−∫𝑇0𝑡𝑟𝑠𝑑𝑠}𝑃^𝛿_0,1𝑇₀,𝑟,𝜆,𝑧;𝑇ℎ^′𝑃₀𝑇₀,𝑟,𝜆,𝑧;𝑇−^𝑇0_𝑡𝑒^{−∫𝑢𝑡𝑟𝑠𝑑𝑠}ℋ₂𝑄₀𝑢,𝑟,𝜆,𝑧;𝑇₀,𝑇𝑑𝑢∣𝑟_𝑡𝜆=𝑟,_𝑡=𝜆,𝑍_𝑡=𝑧. (47)

Synthesizing the above results, we obtain an asymptotic pricing formula for a European option written on the defaultable bond as follows.

Theorem 5. The option price 𝑄(𝑡,𝑟,𝜆,𝑦,𝑧;𝑇₀,𝑇) defined by (34)-(35) is given by 𝑄𝑡,𝑟,𝜆,𝑦,𝑧;𝑇₀,𝑇≈𝑄₀𝑡,𝑟,𝜆,𝑧;𝑇₀+𝑄,𝑇_1,0𝑡,𝑟,𝜆,𝑧;𝑇₀+𝑄,𝑇_0,1𝑡,𝑟,𝜆,𝑧;𝑇₀,𝑇, (48) where 𝑄₀, 𝑄_1,0, and 𝑄_0,1 are given by (41), (44), and (47), respectively.

6. Conclusion

Based upon a reduced form framework of credit risk, this paper investigates the term structure of interest rate derivatives by utilizing an asymptotic expansion method. Mainly, it is focused on the multiscale stochastic volatility of default intensity of the CIR type. Firstly, we obtain an approximation to the value of the defaultable bond as an extension of the known affine solution for the constant volatility. A small correction to the case of the constant volatility has a useful feature that the original model parameters and functions are replaced by the four group parameters. They melt into the group parameters in the averaged form by the ergodic property of the fast mean-reverting OU process, which is a desirable outcome for the purpose of calibration in practical applications. Secondly, we find out the dependence structure of the two scale factors of the stochastic volatility on the bond price and subsequent yield curve, which suggests some flexibility of the multifactor intensity model. This paper, however, has not tested the model to prove this flexibility from empirical evaluation. Furthermore, it is necessary to illustrate the performance of the model setup for the fitting of yield curves. The model calibration and empirical evidence are left to future extensions which also include the pricing of the credit default swap (CDS) spread as opposed to studies based on a Cox process such as [19]. Furthermore, it should be interesting to apply the framework of this paper to multiname intensity models for the pricing of collateralized debt obligation (CDO).

Acknowledgment

The authors thank the anonymous referees for helpful remarks to improve this paper. This study was supported by the National Research Foundation of Korea NRF-2010-0008717.