# Multiscale Analysis on the Pricing of Intensity-Based Defaultable Bonds

Authors: Sun-Hwa Cho; Jeong-Hoon Kim; Yong-Ki Ma
Journal: Journal of Applied Mathematics
 (2013)

---

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


---

## Body

## 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 Section2, 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 ratert, intensity λt, two small positive parameters ϵ and δ, and two processes Yt and Zt representing a fast scale factor and a slow scale factor of the volatility of the intensity λt, respectively, we assume that dynamics of the joint process (rt,λt,Yt,Zt) are given by the following stochastic differential equations (SDEs):
(1)drt=a(r*-rt)dt+σrtdWtr,dλt=a^(λ*-λt)dt+f(Yt,Zt)λtdWtλ,dYt=1ϵ(m-Yt)dt+ν2ϵdWty,dZt=δg(Zt)dt+δh(Zt)dWtz
under a risk-neutral probability measure (or equivalent martingale measure), where Wtr, Wtλ, Wty, and Wtz are standard Brownian motions with a correlation structure given by the coefficients ρrλ, ρry, ρrz, ρλy, ρλz, and ρyz with ρrλ=0. The function f:ℝ2→ℝ+ is assumed to be bounded, smooth, and strictly positive and the functions g and h satisfy the Lipschitz and growth conditions so that the corresponding SDE admits a unique strong solution. We note particularly that the process Yt is an ergodic process whose invariant distribution is given by the Gaussian probability distribution function as
(2)ϕ(y)=12πν2exp(-(y-m)22ν2),y∈ℝ,
which provides an important averaging tool for the unobserved process Yt 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)Pϵ,δ(t,r,λ,y,z;T)=E*[e-∫tT‍(rs+λs)ds∣rt=r,λt=λ,Yt=y,Zt=z]
under the risk-neutral probability measure. Here, the Markov property of the joint process (rt,λt,Yt,Zt) 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
(4)ℒϵ,δPϵ,δ(t,r,λ,y,z;T)=0,t<T,ℒϵ,δ∶=1ϵℒ0+1ϵℒ1+ℒ2+δ𝒦1+δ𝒦2+δϵ𝒦3,Pϵ,δ(t,r,λ,y,z;T)|t=T=1,
where
(5)ℒ0=(m-y)∂∂y+ν2∂2∂y2,ℒ1=2σνρryr∂2∂r∂y+2f(y,z)νρλyλ∂2∂λ∂y,ℒ2=∂∂t+a(r*-r)∂∂r+12σ2r∂2∂r2+a^(λ*-λ)∂∂λ+12f2(y,z)λ∂2∂λ2-(r+λ)I,𝒦1=σh(z)ρrzr∂2∂r∂z+f(y,z)h(z)ρλzλ∂2∂λ∂z,𝒦2=g(z)∂∂z+12h2(z)∂2∂z2,𝒦3=2νh(z)ρyz∂2∂y∂z.
Note that (1/ϵ)ℒ0 is the infinitesimal generator of the OU process Yt. The operator ℒ1 contains the mixed partial derivative due to the correlation between rt and Yt as well as between λt and Yt. ℒ2 is the operator of the canonical two-factor CIR model with volatility at the volatility level f(y,z). 𝒦1 includes the mixed partial derivatives due to the correlation between rt and Zt and between λt and Zt. 𝒦2 is the infinitesimal generator of the process Zt. 𝒦3 holds the mixed partial derivative due to the correlation between Yt and Zt.


## 3. Multiscale Analysis

Since it is difficult to solve the PDE problem (4) itself, we are interested in the following asymptotic expansions:
(6)Pϵ,δ(t,r,λ,y,z;T)=∑j=0∞‍δj/2Pjϵ(t,r,λ,y,z;T),(7)Pjϵ(t,r,λ,y,z;T)=∑i=0∞‍ϵi/2Pi,j(t,r,λ,y,z;T),
so that Pϵ,δ is a series of the general term ϵi/2δj/2Pi,j. Plugging the expansion (6) into (4) leads to P0ϵ and P1ϵ given by the solutions of the PDEs
(8)(1ϵℒ0+1ϵℒ1+ℒ2)P0ϵ(t,r,λ,y,z;T)=0,t<T,P0ϵ(t,r,λ,y,z;T)|t=T=1,(9)(1ϵℒ0+1ϵℒ1+ℒ2)P1ϵ(t,r,λ,y,z;T)=-(𝒦1+1ϵ𝒦3)P0ϵ,t<T,P1ϵ(t,r,λ,y,z;T)|t=T=0,
respectively. From now on, we employ an analytic technique of [13] to approximate the two terms P0ϵ and P1ϵ based on the ergodic property of the Ornstein-Uhlenbeck (OU) process Yt.


### 3.1. The Leading Order Term

Applying the expansion (7) with j=0 to (8) leads to
(10)1ϵℒ0P0,0+1ϵ(ℒ0P1,0+ℒ1P0,0)+(ℒ0P2,0+ℒ1P1,0+ℒ2P0,0)+ϵ(ℒ0P3,0+ℒ1P2,0+ℒ2P1,0)+⋯=0
from which we obtain an affine form of the leading order term P0,0 as follows.Theorem 1.
Assume that the partial derivative ofPi,j with respect to y does not grow as much as ∂Pi,j/∂y~ey2/2 as y goes to infinity. Then the leading order term P0∶=P0,0 of the expansion (7) with j=0 is independent of the fast scale variable y and further it has the affine representation
(11)P0(T-τ,r,λ,z;T)=C˘(τ,z)e-rA˘(τ)-λB˘(τ)
with τ=T-t, A˘(0)=B˘(0)=0 and C˘(0,z)=1, where A˘(τ), B˘(τ), and C˘(τ,z) are given by,
(12)A˘(τ)=2(eγ1τ-1)2γ1+(γ1+a)(eγ1τ-1),B˘(τ)=2(eγ2τ-1)2γ2+(γ2+a^)(eγ2τ-1),C˘(τ,z)=[2γ1e(1/2)(γ1+a)τ2γ1+(γ1+a)(eγ1τ-1)]2ar*/σ2×[2γ2e(1/2)(γ2+a^)τ2γ2+(γ2+a^)(eγ2τ-1)]2a^λ*/σ˘2(z),γ1∶=a2+2σ2,γ2∶=a^2+2σ˘2(z),
respectively. Here, σ˘ is defined by
(13)σ˘2(z)∶=〈f2〉(z)=∫ℝ‍f2(y,z)ϕ(y)dy
in terms of ϕ (the invariant distribution of Yt).Proof.
Multiplying (10) by ϵ and letting ϵ go to zero, we obtain the ordinary differential equation (ODE) ℒ0P0=0. Recalling that the operator ℒ0 is the generator of the OU process Yt, the solution P0 of this ODE must be independent of the y variable due to the assumed growth condition; P0=P0(t,r,λ,z;T). From the O(1/ϵ) terms of (10), we have the PDE ℒ0P1,0+ℒ1P0=0. If we apply the y-independence of P0 to this PDE, then it reduces to the ODE ℒ0P1,0=0 and so P1,0 is independent of y by the same reason as in the case of P0; P1,0=P1,0(t,r,λ,z;T). So far, the two terms P0 and P1,0 do not depend on the current level y of the fast scale volatility driving the process Yt.
One can continue to eliminate the terms of order1, ϵ, ϵ, and so forth. From the O(1) terms of (10), we get ℒ0P2,0+ℒ1P1,0+ℒ2P0=0. This PDE becomes
(14)ℒ0P2,0+ℒ2P0=0
due to the y-independence of P1,0 obtained above. Equation (14) is a Poisson equation for P2,0 with respect to the operator ℒ0 in the y variable. It is well known (cf. [17]) from the Fredholm alternative (solvability condition) that it has a solution only if ℒ2P0 is centered with respect to the invariant distribution of Yt; that is,
(15)〈ℒ2〉P0=0,〈ℒ2〉∶=∂∂t+a(r*-r)∂∂r+12σ2r∂2∂r2+a^(λ*-λ)∂∂λ+12σ˘2(z)λ∂2∂λ2-(r+λ)I.
Noting that (15) is nothing but a PDE corresponding to the intensity with constant volatility replaced by the function σ˘(z), 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 termsP1,0 and P0,1 from the leading order solution P0 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 ofPi,j with respect to y does not grow as much as ∂Pi,j/∂y~ey2/2 as y goes to infinity. The correction term P1,0 is independent of the variable y and P~1,0ϵ∶=ϵP1,0(t,r,λ,z;T) is given by
(16)P~1,0ϵ(T-τ,r,λ,z;T)=D˘(τ,r,λ,z)P0(T-τ,r,λ,z;T)
with D˘(0,r,λ,z)=0, where D˘(τ,r,λ,z) is given by the solution of the PDE
(17)∂D˘∂τ-{a(r*-r)-σ2rA˘(τ)}∂D˘∂r-12σ2r∂2D˘∂r2-{a^(λ*-λ)-σ˘2(z)λB˘(τ)}∂D˘∂λ-12σ˘2(z)λ∂2D˘∂λ2=l(τ,r,λ,z),l(τ,r,λ,z)∶=q(z)λB˘3(τ)-s(z)rA˘(τ)B˘2(τ),q(z)∶=νϵ2ρλyλ〈fθy〉,s(z)∶=νϵ2ρryσr〈θy〉.
Here, A˘(τ) and B˘(τ) are given by (12) and the function θ:ℝ2→ℝ is defined by the solution of
(18)ℒ0θ(y,z)=f2(y,z)-〈f2〉(z)
and θy denotes the partial derivative with respect to the y variable.Proof.
TheO(ϵ) terms of (10) lead to ℒ0P3,0+ℒ1P2,0+ℒ2P1,0=0 which is a Poisson equation for P3,0 and so the solvability condition leads to the PDE
(19)〈ℒ1P2,0+ℒ2P1,0〉=0.
Meanwhile, from (14) and (15), we have
(20)P2,0=-ℒ0-1(ℒ2-〈ℒ2〉)P0+c(t,r,λ,z)
for some function c(t,r,λ,z) that does not rely on the y variable. Plugging (20) into (19), a PDE for P1,0 is obtained as
(21)〈ℒ2〉P1,0=〈ℒ1ℒ0-1(ℒ2-〈ℒ2〉)〉P0
with the terminal condition P1,0(t,r,λ,z;T)|t=T=0. This implies that P1,0 is y-independent. Since we focus on the first-order correction terms to P0, we reset (21) with respect to P~1,0ϵ leading to
(22)〈ℒ2〉P~1,0ϵ=ℋ1P0,ℋ1∶=ϵ〈ℒ1ℒ0-1(ℒ2-〈ℒ2〉)〉.
To calculate the operator ℋ1, we first find
(23)ℒ2-〈ℒ2〉=12(f2(y,z)-〈f2〉(z))∂2∂λ2
from (5) and so the operator ℋ1 is expressed as ℋ1=q(z)λ∂λλλ3+s(z)r∂rλλ3, where q(z) and s(z) are the ones given in the theorem. Then, from (15), Theorem 1 and the change of variables τ=T-t, the PDE (22) leads to
(24)∂P~1,0ϵ∂τ=a(r*-r)∂P~1,0ϵ∂r+12σ2r∂2P~1,0ϵ∂r2+a^(λ*-λ)∂P~1,0ϵ∂λ+12σ˘2(z)λ∂2P~1,0ϵ∂λ2-(r+λ)P~1,0ϵ+C˘(τ,z)e-rA˘(τ)-λB˘(τ)×(q(z)λB˘3(τ)+s(z)rA˘(τ)B˘2(τ))
with the initial condition P~1,0ϵ(T-τ,r,λ,z;T)|τ=0=0. Finally, plugging P~1,0ϵ of the form (16) into (24), we obtain the result of Theorem 2 by direct computation.Similarly, one can derive the correction termP0,1 also as follows.Theorem 3.
Assume that the partial derivative ofPi,j with respect to y does not grow as much as ∂Pi,j/∂y~ey2/2 as y goes to infinity. The first correction term P0,1 does not depend on the variable y and P~0,1δ∶=δP0,1(t,r,λ,z;T) is given by
(25)P~0,1δ(T-τ,r,λ,z;T)=E˘(τ,r,λ,z)P0(T-τ,r,λ,z;T)
with E˘(0,r,λ,z)=0, where E˘(τ,r,λ,z) is given by the solution of the PDE
(26)∂E˘∂τ-{a(r*-r)-σ2rA˘(τ)}∂E˘∂r-12σ2r∂2E˘∂r2-{a^(λ*-λ)-σ˘2(z)λB˘(τ)}∂E˘∂λ-12σ˘2(z)λ∂2E˘∂λ2=m(τ,r,λ,z),m(τ,r,λ,z)∶=1C˘(u(z)rA˘(τ)∂C˘∂z+v(z)λB˘(τ)∂C˘∂z),u(z)∶=-δρrzσh(z),v(z)∶=-δρλzσ-(z)h(z).
Here, A˘(τ), B˘(τ), and C˘(τ,z) are given by (12) and σ- is defined by
(27)σ-(z)∶=〈f〉(z)=∫ℝ‍f(y,z)ϕ(y)dy
in terms of ϕ (the invariant distribution of Yt).Proof.
Applying the expansion (7) with j=0 and j=1 to (9) results in
(28)1ϵℒ0P0,1+1ϵ(ℒ0P1,1+ℒ1P0,1)+(ℒ0P2,1+ℒ1P1,1+ℒ2P0,1)+ϵ(ℒ0P3,1+ℒ1P2,1+ℒ2P1,1)+⋯=-1ϵ𝒦3P0-(𝒦1P0+𝒦3P1,0)-ϵ(𝒦1P1,0+𝒦3P2,0)+⋯.
Then, from the O(1/ϵ) term of (28), P0,1 is y-independent as the solution of the ODE ℒ0P0,1=0 under the assumed growth condition. Since both P0 (obtained in Theorem 1) and P0,1 are y-independent, from the O(1/ϵ) terms of (28), we have the ODE ℒ0P1,1=0. So, by the same reason as in the case of P0,1=0, P1,1 is also independent of y. Hence, the two terms P0,1 and P1,1 do not depend on the current level y of the fast scale volatility of intensity.
One can continue to remove the terms of order1, ϵ, ϵ, and so forth. From the O(1) terms and the y-independence of P1,0 and P1,1, we have the PDE ℒ0P2,1+ℒ2P0,1+𝒦1P0=0. This is a Poisson equation for P2,1 with respect to the operator ℒ0 in the y variable and so it has a solution only if 〈ℒ2〉P0,1=-〈𝒦1〉P0 holds. If we reset this PDE with respect to P~0,1δ, then we have
(29)〈ℒ2〉P~0,1δ=ℋ2P0,ℋ2∶=-δ〈𝒦1〉,
where the operator ℋ2 is the same as ℋ2=u(z)r∂rz2+v(z)λ∂λz2 from the definition of 𝒦1 given by (5).
Using the operator (15), Theorem 1, and the change of variable τ=T-t, we write the PDE (29) as follows:
(30)∂P~0,1δ∂τ=a(r*-r)∂P~0,1δ∂r+12σ2r∂2P~0,1δ∂r2+a^(λ*-λ)∂P~0,1δ∂λ+12σ˘2(z)λ∂2P~0,1δ∂λ2-(r+λ)P~0,1δ+e-rA˘(τ)-λB˘(τ)(u(z)rA˘(τ)∂C˘∂z+e-rA˘(τ)-λB˘(τ)+v(z)λB˘(τ)∂C˘∂z),
where the initial condition P~0,1δ(T-τ,r,λ,z;T)|τ=0=0 is satisfied. Plugging the form (25) into (30), we obtain the result of Theorem 3 by direct computation.


### 3.3. Accuracy

Synthesizing Theorems1, 2, and 3, we obtain the following asymptotic representation of the price of the defaultable bond at time t:
(31)Pϵ,δ(t,r,λ,y,z;T)≈P~ϵ,δ∶=P0+P~1,0ϵ+P~0,1δ=(1+D˘(τ,r,λ,z)+E˘(τ,r,λ,z))×C˘(τ,z)e-rA˘(τ)-λB˘(τ),
where A˘(τ), B˘(τ), C˘(τ,z), D˘(τ,r,λ,z), and E˘(τ,r,λ,z) are given by (12), (17), and (26), respectively. We note that the function g(z) (the drift term of the slow scale variation) in (1) does not affect the present form of the approximated bond price P~ϵ,δ due to the order of δ in front of g(z) in (1). It should influence the next order approximation.In this approximation, the leading order price is determined by the functionsA˘(τ), B˘(τ), and C˘(τ,z) and the first-order corrections are given by D˘(τ,r,λ,z) and E˘(τ,r,λ,z). Here, D˘(τ,r,λ,z) is determined by the two group parameters q(z) and s(z) and E˘(τ,r,λ,z) is determined by u(z) and v(z). These four group parameters absorb the original parameters and functions related to Yt and Zt. The group parameter q(z) absorbs ϵ, ν, ρλy, and f while s(z) absorbs ϵ, ν, ρry, and f. The group parameter u(z) absorbs δ, ρrz, and h while v(z) absorbs δ, ρλz, and h. Particularly, the functions f, g and h 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 Yt 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.
LetP~ϵ,δ be defined by (31). Then for any fixed (t,r,λ,y,z) there exists a positive constant C, which depends on (t,r,λ,y,z) but not on ϵ and δ, such that
(32)Pϵ,δ(t,r,λ,y,z;T)-P~ϵ,δ(t,r,λ,z;T)≤C(ϵ+δ)
holds for all 0<ϵ,  δ≤1.Proof.
To prove the inequality (32), we first define a residual Rϵ,δ by writing Pϵ,δ as
(33)Pϵ,δ=∑i,j=0i,j=3‍ϵi/2δj/2Pi,j-Rϵ,δ.
The second step is to compute ℒϵ,δRϵ,δ using the properties of Pϵ,δ,  P0,…,P3,3 obtained above. The last step is to write Rϵ,δ 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.


## 3.1. The Leading Order Term

Applying the expansion (7) with j=0 to (8) leads to
(10)1ϵℒ0P0,0+1ϵ(ℒ0P1,0+ℒ1P0,0)+(ℒ0P2,0+ℒ1P1,0+ℒ2P0,0)+ϵ(ℒ0P3,0+ℒ1P2,0+ℒ2P1,0)+⋯=0
from which we obtain an affine form of the leading order term P0,0 as follows.Theorem 1.
Assume that the partial derivative ofPi,j with respect to y does not grow as much as ∂Pi,j/∂y~ey2/2 as y goes to infinity. Then the leading order term P0∶=P0,0 of the expansion (7) with j=0 is independent of the fast scale variable y and further it has the affine representation
(11)P0(T-τ,r,λ,z;T)=C˘(τ,z)e-rA˘(τ)-λB˘(τ)
with τ=T-t, A˘(0)=B˘(0)=0 and C˘(0,z)=1, where A˘(τ), B˘(τ), and C˘(τ,z) are given by,
(12)A˘(τ)=2(eγ1τ-1)2γ1+(γ1+a)(eγ1τ-1),B˘(τ)=2(eγ2τ-1)2γ2+(γ2+a^)(eγ2τ-1),C˘(τ,z)=[2γ1e(1/2)(γ1+a)τ2γ1+(γ1+a)(eγ1τ-1)]2ar*/σ2×[2γ2e(1/2)(γ2+a^)τ2γ2+(γ2+a^)(eγ2τ-1)]2a^λ*/σ˘2(z),γ1∶=a2+2σ2,γ2∶=a^2+2σ˘2(z),
respectively. Here, σ˘ is defined by
(13)σ˘2(z)∶=〈f2〉(z)=∫ℝ‍f2(y,z)ϕ(y)dy
in terms of ϕ (the invariant distribution of Yt).Proof.
Multiplying (10) by ϵ and letting ϵ go to zero, we obtain the ordinary differential equation (ODE) ℒ0P0=0. Recalling that the operator ℒ0 is the generator of the OU process Yt, the solution P0 of this ODE must be independent of the y variable due to the assumed growth condition; P0=P0(t,r,λ,z;T). From the O(1/ϵ) terms of (10), we have the PDE ℒ0P1,0+ℒ1P0=0. If we apply the y-independence of P0 to this PDE, then it reduces to the ODE ℒ0P1,0=0 and so P1,0 is independent of y by the same reason as in the case of P0; P1,0=P1,0(t,r,λ,z;T). So far, the two terms P0 and P1,0 do not depend on the current level y of the fast scale volatility driving the process Yt.
One can continue to eliminate the terms of order1, ϵ, ϵ, and so forth. From the O(1) terms of (10), we get ℒ0P2,0+ℒ1P1,0+ℒ2P0=0. This PDE becomes
(14)ℒ0P2,0+ℒ2P0=0
due to the y-independence of P1,0 obtained above. Equation (14) is a Poisson equation for P2,0 with respect to the operator ℒ0 in the y variable. It is well known (cf. [17]) from the Fredholm alternative (solvability condition) that it has a solution only if ℒ2P0 is centered with respect to the invariant distribution of Yt; that is,
(15)〈ℒ2〉P0=0,〈ℒ2〉∶=∂∂t+a(r*-r)∂∂r+12σ2r∂2∂r2+a^(λ*-λ)∂∂λ+12σ˘2(z)λ∂2∂λ2-(r+λ)I.
Noting that (15) is nothing but a PDE corresponding to the intensity with constant volatility replaced by the function σ˘(z), 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 termsP1,0 and P0,1 from the leading order solution P0 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 ofPi,j with respect to y does not grow as much as ∂Pi,j/∂y~ey2/2 as y goes to infinity. The correction term P1,0 is independent of the variable y and P~1,0ϵ∶=ϵP1,0(t,r,λ,z;T) is given by
(16)P~1,0ϵ(T-τ,r,λ,z;T)=D˘(τ,r,λ,z)P0(T-τ,r,λ,z;T)
with D˘(0,r,λ,z)=0, where D˘(τ,r,λ,z) is given by the solution of the PDE
(17)∂D˘∂τ-{a(r*-r)-σ2rA˘(τ)}∂D˘∂r-12σ2r∂2D˘∂r2-{a^(λ*-λ)-σ˘2(z)λB˘(τ)}∂D˘∂λ-12σ˘2(z)λ∂2D˘∂λ2=l(τ,r,λ,z),l(τ,r,λ,z)∶=q(z)λB˘3(τ)-s(z)rA˘(τ)B˘2(τ),q(z)∶=νϵ2ρλyλ〈fθy〉,s(z)∶=νϵ2ρryσr〈θy〉.
Here, A˘(τ) and B˘(τ) are given by (12) and the function θ:ℝ2→ℝ is defined by the solution of
(18)ℒ0θ(y,z)=f2(y,z)-〈f2〉(z)
and θy denotes the partial derivative with respect to the y variable.Proof.
TheO(ϵ) terms of (10) lead to ℒ0P3,0+ℒ1P2,0+ℒ2P1,0=0 which is a Poisson equation for P3,0 and so the solvability condition leads to the PDE
(19)〈ℒ1P2,0+ℒ2P1,0〉=0.
Meanwhile, from (14) and (15), we have
(20)P2,0=-ℒ0-1(ℒ2-〈ℒ2〉)P0+c(t,r,λ,z)
for some function c(t,r,λ,z) that does not rely on the y variable. Plugging (20) into (19), a PDE for P1,0 is obtained as
(21)〈ℒ2〉P1,0=〈ℒ1ℒ0-1(ℒ2-〈ℒ2〉)〉P0
with the terminal condition P1,0(t,r,λ,z;T)|t=T=0. This implies that P1,0 is y-independent. Since we focus on the first-order correction terms to P0, we reset (21) with respect to P~1,0ϵ leading to
(22)〈ℒ2〉P~1,0ϵ=ℋ1P0,ℋ1∶=ϵ〈ℒ1ℒ0-1(ℒ2-〈ℒ2〉)〉.
To calculate the operator ℋ1, we first find
(23)ℒ2-〈ℒ2〉=12(f2(y,z)-〈f2〉(z))∂2∂λ2
from (5) and so the operator ℋ1 is expressed as ℋ1=q(z)λ∂λλλ3+s(z)r∂rλλ3, where q(z) and s(z) are the ones given in the theorem. Then, from (15), Theorem 1 and the change of variables τ=T-t, the PDE (22) leads to
(24)∂P~1,0ϵ∂τ=a(r*-r)∂P~1,0ϵ∂r+12σ2r∂2P~1,0ϵ∂r2+a^(λ*-λ)∂P~1,0ϵ∂λ+12σ˘2(z)λ∂2P~1,0ϵ∂λ2-(r+λ)P~1,0ϵ+C˘(τ,z)e-rA˘(τ)-λB˘(τ)×(q(z)λB˘3(τ)+s(z)rA˘(τ)B˘2(τ))
with the initial condition P~1,0ϵ(T-τ,r,λ,z;T)|τ=0=0. Finally, plugging P~1,0ϵ of the form (16) into (24), we obtain the result of Theorem 2 by direct computation.Similarly, one can derive the correction termP0,1 also as follows.Theorem 3.
Assume that the partial derivative ofPi,j with respect to y does not grow as much as ∂Pi,j/∂y~ey2/2 as y goes to infinity. The first correction term P0,1 does not depend on the variable y and P~0,1δ∶=δP0,1(t,r,λ,z;T) is given by
(25)P~0,1δ(T-τ,r,λ,z;T)=E˘(τ,r,λ,z)P0(T-τ,r,λ,z;T)
with E˘(0,r,λ,z)=0, where E˘(τ,r,λ,z) is given by the solution of the PDE
(26)∂E˘∂τ-{a(r*-r)-σ2rA˘(τ)}∂E˘∂r-12σ2r∂2E˘∂r2-{a^(λ*-λ)-σ˘2(z)λB˘(τ)}∂E˘∂λ-12σ˘2(z)λ∂2E˘∂λ2=m(τ,r,λ,z),m(τ,r,λ,z)∶=1C˘(u(z)rA˘(τ)∂C˘∂z+v(z)λB˘(τ)∂C˘∂z),u(z)∶=-δρrzσh(z),v(z)∶=-δρλzσ-(z)h(z).
Here, A˘(τ), B˘(τ), and C˘(τ,z) are given by (12) and σ- is defined by
(27)σ-(z)∶=〈f〉(z)=∫ℝ‍f(y,z)ϕ(y)dy
in terms of ϕ (the invariant distribution of Yt).Proof.
Applying the expansion (7) with j=0 and j=1 to (9) results in
(28)1ϵℒ0P0,1+1ϵ(ℒ0P1,1+ℒ1P0,1)+(ℒ0P2,1+ℒ1P1,1+ℒ2P0,1)+ϵ(ℒ0P3,1+ℒ1P2,1+ℒ2P1,1)+⋯=-1ϵ𝒦3P0-(𝒦1P0+𝒦3P1,0)-ϵ(𝒦1P1,0+𝒦3P2,0)+⋯.
Then, from the O(1/ϵ) term of (28), P0,1 is y-independent as the solution of the ODE ℒ0P0,1=0 under the assumed growth condition. Since both P0 (obtained in Theorem 1) and P0,1 are y-independent, from the O(1/ϵ) terms of (28), we have the ODE ℒ0P1,1=0. So, by the same reason as in the case of P0,1=0, P1,1 is also independent of y. Hence, the two terms P0,1 and P1,1 do not depend on the current level y of the fast scale volatility of intensity.
One can continue to remove the terms of order1, ϵ, ϵ, and so forth. From the O(1) terms and the y-independence of P1,0 and P1,1, we have the PDE ℒ0P2,1+ℒ2P0,1+𝒦1P0=0. This is a Poisson equation for P2,1 with respect to the operator ℒ0 in the y variable and so it has a solution only if 〈ℒ2〉P0,1=-〈𝒦1〉P0 holds. If we reset this PDE with respect to P~0,1δ, then we have
(29)〈ℒ2〉P~0,1δ=ℋ2P0,ℋ2∶=-δ〈𝒦1〉,
where the operator ℋ2 is the same as ℋ2=u(z)r∂rz2+v(z)λ∂λz2 from the definition of 𝒦1 given by (5).
Using the operator (15), Theorem 1, and the change of variable τ=T-t, we write the PDE (29) as follows:
(30)∂P~0,1δ∂τ=a(r*-r)∂P~0,1δ∂r+12σ2r∂2P~0,1δ∂r2+a^(λ*-λ)∂P~0,1δ∂λ+12σ˘2(z)λ∂2P~0,1δ∂λ2-(r+λ)P~0,1δ+e-rA˘(τ)-λB˘(τ)(u(z)rA˘(τ)∂C˘∂z+e-rA˘(τ)-λB˘(τ)+v(z)λB˘(τ)∂C˘∂z),
where the initial condition P~0,1δ(T-τ,r,λ,z;T)|τ=0=0 is satisfied. Plugging the form (25) into (30), we obtain the result of Theorem 3 by direct computation.


## 3.3. Accuracy

Synthesizing Theorems1, 2, and 3, we obtain the following asymptotic representation of the price of the defaultable bond at time t:
(31)Pϵ,δ(t,r,λ,y,z;T)≈P~ϵ,δ∶=P0+P~1,0ϵ+P~0,1δ=(1+D˘(τ,r,λ,z)+E˘(τ,r,λ,z))×C˘(τ,z)e-rA˘(τ)-λB˘(τ),
where A˘(τ), B˘(τ), C˘(τ,z), D˘(τ,r,λ,z), and E˘(τ,r,λ,z) are given by (12), (17), and (26), respectively. We note that the function g(z) (the drift term of the slow scale variation) in (1) does not affect the present form of the approximated bond price P~ϵ,δ due to the order of δ in front of g(z) in (1). It should influence the next order approximation.In this approximation, the leading order price is determined by the functionsA˘(τ), B˘(τ), and C˘(τ,z) and the first-order corrections are given by D˘(τ,r,λ,z) and E˘(τ,r,λ,z). Here, D˘(τ,r,λ,z) is determined by the two group parameters q(z) and s(z) and E˘(τ,r,λ,z) is determined by u(z) and v(z). These four group parameters absorb the original parameters and functions related to Yt and Zt. The group parameter q(z) absorbs ϵ, ν, ρλy, and f while s(z) absorbs ϵ, ν, ρry, and f. The group parameter u(z) absorbs δ, ρrz, and h while v(z) absorbs δ, ρλz, and h. Particularly, the functions f, g and h 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 Yt 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.
LetP~ϵ,δ be defined by (31). Then for any fixed (t,r,λ,y,z) there exists a positive constant C, which depends on (t,r,λ,y,z) but not on ϵ and δ, such that
(32)Pϵ,δ(t,r,λ,y,z;T)-P~ϵ,δ(t,r,λ,z;T)≤C(ϵ+δ)
holds for all 0<ϵ,  δ≤1.Proof.
To prove the inequality (32), we first define a residual Rϵ,δ by writing Pϵ,δ as
(33)Pϵ,δ=∑i,j=0i,j=3‍ϵi/2δj/2Pi,j-Rϵ,δ.
The second step is to compute ℒϵ,δRϵ,δ using the properties of Pϵ,δ,  P0,…,P3,3 obtained above. The last step is to write Rϵ,δ 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 factorD˘ and the slow term factor E˘ 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.Table1 shows the average cumulative default rates of bond issuers for the period of 1983–2009. It says, for instance, that over a ten-year B-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 t from the cumulative probability of default by time t. For example, the average 20-year default intensity for B-rated issuers is given by -(1/20)ln(1-0.65493)=0.0532.Table 1
Average cumulative default rates (%), 1983–2009 (source: Moody's).
Rating
Year 1
Year 2
Year 5
Year 7
Year 10
Year 15
Year 20
Aaa
0
0.016
0.086
0.182
0.187
0.187
0.187
Aa
0.024
0.066
0.247
0.318
0.408
0.989
2.583
A
0.057
0.187
0.806
1.297
2.099
3.521
6.536
Baa
0.196
0.543
2.027
3.130
4.815
8.719
12.603
Ba
1.209
3.434
11.444
15.600
21.128
30.467
38.152
B
4.550
10.519
26.24
35.054
44.982
57.136
65.493
Caa
15.383
26.969
50.339
57.783
71.993
82.434
86.669
Ca-C
36.207
48.440
70.176
74.757
78.014
78.014
—
All rated
1.761
3.620
8.046
9.991
12.100
14.990
18.002Based on Table1 and the data [13], we choose the parameter values a=1.2,  a^=1.3,  σ=0.06,  σ˘2=0.03,  r*=0.1,  λ=0.0455,  and  λ*=0.0532. Maturity runs from 0 to 20 years and the interest rate is fixed as r=0.07. Figures 1 and 2 show the bond price correction term structure depending on the group parameters (q,s) and (u,v), respectively, for the fast term factor D˘ and the slow term factor E˘. 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 factorD˘ is shown by dotted, dashed, and sold lines as functions of time to maturity, respectively, for (q,s)=(0.2,0.1), (0.3,0.4), and (0.6,0.7).Figure 2
The price correction term structure factorE˘ is shown by dotted, dashed and sold lines as functions of time to maturity, respectively, for (u,v)=(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 notationsT and T0,  T0<T, 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 T0. The short rate process rt, the intensity process λt, the fast volatility process Yt of the intensity, and the slow volatility process Zt of the intensity are given by the SDEs (1).The option price at timet for an observed short rate rt=r, an intensity level λt=λ, a fast volatility level Yt=y, and a slow volatility level Zt=z, denoted by Q(t,r,λ,y,z;T0,T), is defined by
(34)Q(t,r,λ,y,z;T0,T)=E*[e-∫tT0‍rsdsh(ℬ(T0,T))∣rt=r,λt=λ,Yt=y,Zt=z]
under an equivalent martingale measure, where the bond price ℬ(T0,T) is given by
(35)ℬ(T0,T)=Pϵ,δ(T0,rT0,λT0,YT0,ZT0;T)=E*[e-∫T0T‍(rs+λs)ds∣rT0,λT0,YT0,ZT0]
and h is the payoff function of the option at time T0. Here, it is assumed that the payoff function h 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 h is not differentiable at the exercise price in the case of classical European call or put option. The smoothness assumption on h, however, can be removed similarly to the argument in [18] or [12]. So, for simplicity, h is assumed to be a smooth function.By the four-dimensional Feynman-Kac formula, we transform the above integral form (34) into the PDE
(36)∂Q∂t+a(r*-r)∂Q∂r+12σ2r∂2Q∂r2+a^(λ*-λ)∂Q∂λ+12f2(y,z)λ∂2Q∂λ2+1ϵ(m-y)∂Q∂y+ν2ϵ∂2Q∂y2+δg(z)∂Q∂z+12δh2(z)∂2Q∂z2+ρryσrν2ϵ∂2Q∂r∂y+ρrzσrδh(z)∂2Q∂r∂z+ρλyf(y,z)λν2ϵ∂2Q∂λ∂y+ρλzf(y,z)λδh(z)∂2Q∂λ∂z+ρyzν2ϵδh(z)∂2Q∂y∂z-(r+λ)Q=0
with the terminal condition Q(t,r,λ,y,z;T0,T)|t=T0=h(Pϵ,δ(T0,r,λ,y,z;T)). Then, keeping the notation used in Section 3 but with a different terminal condition, the option price Q(t,r,λ,y,z;T0,T) is given by the solution of the PDE problem
(37)ℒϵ,δQ(t,r,λ,y,z;T0,T)=0,Q(t,r,λ,y,z;T0,T)|t=T0=h(Pϵ,δ(T0,r,λ,y,z;T)).
Taking the asymptotic expansions
(38)Q=Q0ϵ+δQ1ϵ+δQ2ϵ+⋯,(39)Qkϵ=Q0,k+ϵQ1,k+ϵQ2,k+ϵ3/2Q3,k+⋯,k=0,1,2,…,(40)Q(t,r,λ,y,z;T0,T)|t=T0=h(P0(T0,r,λ,z;T))+P~1,0ϵ(T0,r,λ,z;T)h′(P0(T0,r,λ,z;T))+P~0,1δ(T0,r,λ,z;T)h′(P0(T0,r,λ,z;T))+⋯,
where P0, P~1,0ϵ, and P~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 order1/ϵ and 1/ϵ in the asymptotic PDE (37) yield the y-independence of Q0,0 (in brief, Q0), Q1,0, and Q0,1 under the same growth condition as in Section 3.TheO(1) terms in (37) give a Poisson equation for Q2,0 from which the solvability condition 〈ℒ2〉Q0=0 has to be satisfied. From (40), the corresponding terminal condition is given by Q0|t=T0=h(P0(T0,r,λ,z;T)). Then we have
(41)Q0(t,r,λ,z;T0,T)=E*[e-∫tT0‍rsdsh(P0(T0,r,λ,z;T))∣rt=r,λt=λ,Zt=z],
where P0(T0,r,λ,z;T) is given by (11) in Theorem 1 at time t=T0.TheO(ϵ) terms in (37) lead to a Poisson equation for Q3,0, where the solvability condition is given by 〈ℒ1Q2,0+ℒ2Q1,0〉=0. If we put Q~1,0∶=ϵQ1,0, then this solvability condition leads to
(42)〈ℒ2〉Q~1,0(t,r,λ,z;T0,T)=ℋ1Q0(t,r,λ,z;T0,T)
with a terminal condition given by
(43)Q~1,0(t,r,λ,z;T0,T)|t=T0=P~1,0ϵ(T0,r,λ,z;T)h′(P0(T0,r,λ,z;T))
from (40), where the operator ℋ1 is given by (22) and P~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:
(44)Q~1,0(t,r,λ,z;T0,T)=E*[∫tT0‍e-∫tT0‍rsdsP~1,0ϵ(T0,r,λ,z;T)h′(P0(T0,r,λ,z;T))-∫tT0‍e-∫tu‍rsdsℋ1Q0(u,r,λ,z;T0,T)du∣rt=r,λt=λ,Zt=z∫tT0‍].
Now, the O(δ) terms in (37) lead to a Poisson equation for Q2,1 such that the solvability condition 〈ℒ2Q0,1〉=-〈𝒦1Q0〉 holds. If we let Q~0,1∶=δQ0,1, then this solvability condition leads to
(45)〈ℒ2〉Q~0,1(t,r,λ,z;T0,T)=ℋ2Q0(t,r,λ,z;T0,T),
where the operator ℋ2 is given by (29). Also, from the terminal condition (40), we have
(46)Q~0,1(t,r,λ,z;T0,T)|t=T0=P~0,1δ(T0,r,λ,z;T)h′(P0(T0,r,λ,z;T)),
where P~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:
(47)Q~0,1(t,r,λ,z;T0,T)=E*[∫tT0‍e-∫tT0‍rsdsP~0,1δ(T0,r,λ,z;T)h′(P0(T0,r,λ,z;T))-∫tT0‍e-∫tu‍rsdsℋ2Q0(u,r,λ,z;T0,T)du∣rt=r,λt=λ,Zt=z∫tT0‍].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 priceQ(t,r,λ,y,z;T0,T) defined by (34)-(35) is given by
(48)Q(t,r,λ,y,z;T0,T)≈Q0(t,r,λ,z;T0,T)+Q~1,0(t,r,λ,z;T0,T)+Q~0,1(t,r,λ,z;T0,T),
where Q0, Q~1,0, and Q~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).



---

*Source: 287425-2013-05-16.xml*