European Journal of Operational Research 185 (2008) 1462–1476 www.elsevier.com/locate/ejor
Calibration of the default probability model Alexander Kreinin *, Ahmed Nagi
1
Algorithmics Inc., 185 Spadina Avenue, Toronto, Ont., Canada M5T 2C6 Received 1 November 2003; accepted 1 November 2004 Available online 28 November 2006
Abstract In this paper, we study the calibration problem for the Merton–Vasicek default probability model [Robert Merton, On the pricing of corporate debt: the risk structure of interest rate, Journal of Finance 29 (1974) 449–470]. We derive conditions that guarantee existence and uniqueness of the solution. Using analytical properties of the model, we propose a fast calibration procedure for the conditional default probability model in the integrated market and credit risk framework. Our solution allows one to avoid numerical integration problems as well as problems related to the numerical solution of the nonlinear equations. Ó 2006 Elsevier B.V. All rights reserved. Keywords: Credit risk; Default probability; Merton–Vasicek model; Calibration
1. Introduction Calibration of the default probability model in the integrated market and credit risk framework [7] is one of the problems that require intensive computations. Indeed, this problem should be solved separately for each obligor. The number of obligors contained in a corporate portfolio may be very large. Taking into account that the straightforward implementation of the calibration algorithm (algorithm NINS) contains numerical integration and numerical solution of a nonlinear equation [11], we conclude that an efficient solution to this problem becomes critical for an industrial implementation of this model. An alternative solution can be based on reduction of the calibration problem to computation of the bivariate normal integral [14]. The latter approach is faster than NINS algorithm [11] but also requires numerical solution of the nonlinear equation for the parameters of the model. In this paper, we describe an efficient approach to the calibration of the default probability model based on analytical properties of the calibration equation in the integrated market and credit risk framework. 2 In this paper, we consider the calibration problem for the one step Merton–Vasicek default model. This model is
*
1 2
Corresponding author. E-mail address:
[email protected] (A. Kreinin). Present address: Citigroup, 250 West Street, New York, NY 10013, USA. The integrated market and credit risk model is described in details in [7,6].
0377-2217/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2004.11.029
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
1463
described, in particular, in [9,13,12,4–6,8]. This model is actively used in one-factor setting for estimation of asset correlations [2,3]. The calibration procedure studied in this paper is based on matching the first two moments of the conditional default probabilities. We find an approximate solution to the calibration problem that does not require numerical integration or numerical solution of the nonlinear equation. Our solution provides sufficiently small level of relative error acceptable for practical implementation. Our approach allows us to capture asymptotic behaviour of the parameters when the correlation of the obligor’s creditworthiness index and the credit driver is close to 0 as well as close to ±1. In the latter case application of maximum likelihood method usually leads to numerically unstable results. This paper is organized as follows. In Section 2, the model for conditional default probability of the obligor is introduced. In Section 3.1, we study the moments of the conditional default probabilities. In Section 4, we construct three different approximations of the second moment of the default probability covering different areas of the parameters value. These approximations are obtained using the results of Section 3.1. In Section 5, we propose a calibration procedure based on the approximations studied in Section 4. In Section 6, we analyze the accuracy of our solution and compare the computation time of the straightforward calibration procedure and that proposed in Section 5. 2. The model In the integrated market and credit risk model, the default event is governed by a credit worthiness index, Yt. This index is described by the random process Y t ¼ b nt þ a e t ;
t ¼ 1; 2; . . . ;
a2 þ b2 ¼ 1;
ð1Þ
where both nt and et are independent canonical processes [6], i.e. the random variables nt and et all have standard normal distribution Nð0; 1Þ, t = 1, 2, . . . The processes nt and et are interpreted as follows: nt represents a systematic component of the market risk (including credit drivers) and et represents an idiosyncratic component. The default time, s, in this model is represented as a first hitting time of the process Yt s ¼ inf ft : Y t < bðtÞg; t>0
where the boundary b( Æ ) is either a deterministic continuous function in the continuous-time case or a deterministic sequence in the discrete-time case. Consider the default event (Y1 < b) in a scenario in which the value of the systematic component is equal to x. Then from (1) we have for t = 1, Y1 = b Æ x + a Æ e1, and inR the case of default e1 < bbx . a z u2 =2 ffi 1 Denote p = Pr{Y1 < b}. Then p = U(b), where UðzÞ ¼ p1ffiffiffi e du is the cumulative distribution func2p tion of a standard normal random variable. The conditional probability of default, P ðn1 Þ ¼ PrðY 1 < bjn1 Þ; therefore, is
b b n1 P ðn1 Þ ¼ U : a
ð2Þ
Since n1 has a standard normal distribution, we can consider the distribution of the conditional default probability g = P(n1). We have from (2) b b n1 b a U1 ðvÞ Prfg < vg ¼ Pr U ; where 0 6 v 6 1: ð3Þ < v ¼ Pr n1 > b a Denote the cumulative distribution function of g by F(v). Then from (3) we derive b a U1 ðvÞ F ðvÞ ¼ U ; b where UðÞ ¼ 1 UðÞ.
ð4Þ
1464
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
Let us denote by fb(v) the probability density function of g. Then from (4) we have 1 2 2 2 1 1 fb ðvÞ ¼ k exp ðc 2 ckU ðvÞ þ ðk 1Þ ðU ðvÞÞ Þ ; 2
ð5Þ
where a k¼ ; b
b c¼ : b
It follows from (5) that for sufficiently small unconditional default probabilities p < 0.5 we have c < 0 and then 1; if k 6 1; lim fb ðvÞ ¼ v!0 0; if k > 1 and lim fb ðvÞ ¼ v!1
1;
if k < 1;
0;
if k P 1:
If p = 0.5 then c = 0. In this case we have fb(v) = 1, if k = 1; if k < 1 then lim fb ðvÞ ¼ þ1; v!0
and
lim fb ðvÞ ¼ þ1 v!1
If k > 1 then lim fb ðvÞ ¼ 0;
and
v!0
lim fb ðvÞ ¼ 0: v!1
Fig. 1 illustrates the variety of shapes of the family fb of the default probability density consistent with the relations for lim fb ðvÞ and lim fb ðvÞ. v!0
v!1
p= 0.37 β=0.1
p= 0.37 β = 0.707106
10
2.5
8
2 f (v)
3
6
1.5
β
fβ(v)
12
4
1
2
0.5
0
0
0.2
0.4
0.6
0.8
0
1
0
0.2
0.4
0.6
v
0.8
1
0.8
1
v
p= 0.37 β=0.95
p= 0.05 β=0.707106
5
15
4 10 β
f (v)
fβ(v)
3 2
5 1 0
0
0.2
0.4
0.6 v
0.8
1
0
0
0.2
0.4
0.6 v
Fig. 1. Density of the default probability distribution.
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
1465
3. The calibration problem The model described by Eq. (2) contains two unknown parameters, b and b. Therefore, if we can find two independent equations linking these parameters to some statistical characteristics of g = P(x) we will be able to solve the calibration problem. In this Section, we apply the idea of matching the first two moments of the random variable g. Thus, our analysis will focus on the solution of the calibration equations Egðb; bÞ ¼ p;
ð6Þ
Eg2 ðb; bÞ ¼ Q;
ð7Þ
where Q = p2 + r2 and r is an empirical standard deviation of the conditional default probability. Using analytical properties of the function Qðb; bÞ ¼ Eg2 ðb; bÞ we will prove that the calibration problems (6) and (7) has a unique solution b ¼ bðp; QÞ;
b ¼ bðp; QÞ;
under some natural constraints on p and Q. In Section 3.1, we prove the main analytical result – integral representation of the function Q(b, b). The proof is based on four lemmas that describe the properties of the second moment of the default probability. The derivation of these four statements is placed in Appendix A. 3.1. The second moment of the default probability The first important property of the unconditional default probability is invariance of Eg with respect to the parameter b. Lemma 1. The expectation of the conditional default probability satisfies the relation Egðb; bÞ ¼ p;
ð8Þ
for all b, 1 6 b 6 1, where p = U1(b) is the unconditional default probability. Proof. The proof of (8) is based on the following statement. Proposition 2. Let n be a random variable having a standard normal distribution, n Nð0; 1Þ. Then the expectation of the random variable U(a Æ n + c) equals c EUða n þ cÞ ¼ U pffiffiffiffiffiffiffiffiffiffiffiffiffi : ð9Þ 1 þ a2 Proposition 2 is known. For the sake of completeness we give a short proof of (9) in the Appendix A. Let us now derive (8) from Proposition 2. The expectation of the conditional default probability, Egðb; bÞ, satisfies the relation bbn Egðb; bÞ ¼ EU ; a
a2 þ b2 ¼ 1;
UðbÞ ¼ p;
where n is a random variable having a standard normal distribution. Then, applying (9), we obtain 0 1 b
B C a ffiA ¼ UðbÞ ¼ p; Egðb; bÞ ¼ U@qffiffiffiffiffiffiffiffiffiffiffi b2 1 þ a2 as was to be proved.
1466
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
Lemma 1 implies that the equation Egðb; bÞ ¼ p; has a unique solution for the parameter b, b = U1(p). Therefore, to solve the calibration problem one has to find a solution of the equation Eg2 ðb; bÞ ¼ Q given b, where Q is the empirical second moment of the conditional default probability. The second moment, Qðb; bÞ ¼ Eg2 ðb; bÞ, of the conditional default probability satisfies the relation3 0 1 Z 1 B b bx C Qðb; bÞ ¼ U2 @qffiffiffiffiffiffiffiffiffiffiffiffiffiAdUðxÞ; 1 < b < 1; 1 6 b 6 1: ð10Þ 2 1 1b Let us consider analytical properties of the function Q(b, b). Denote sffiffiffiffiffiffiffiffiffiffiffiffiffi 1 b2 : kðbÞ ¼ 1 þ b2 Theorem 3. The function Q(b, b) has the following integral representation Z kðbÞ exp b2 ðt2 þ1Þ 2 1 dt: Qðb; bÞ ¼ UðbÞ p 0 t2 þ 1 The partial derivative
oQðb;bÞ ob
ð11Þ
equals
2 oQðb; bÞ 1 b b ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffi e 1þb2 : ob p 1 b2 1 þ b2
ð12Þ
Before sketching the proof of Theorem 3, we formulate an important corollary that states monotonicity of the function Q. Corollary 4. For fixed b, the function Q(b, b) is a monotonically increasing function of b for b P 0. The proof of Corollary 4 follows directly from (12). Proof of Theorem 3. The proof is based on the following statements. Lemma 5. The partial derivative of the function Q(b, b) satisfies the relation oQðb; bÞ ¼ 2uðbÞUðb kðbÞÞ; ob pffiffiffiffiffiffi where uðxÞ ¼ expðx2 =2Þ= 2p is the density function of the standard normal distribution. Formula (13) implies that the function Q(b, b) can be represented in the form Z b Qðb; bÞ ¼ Qð0; bÞ þ 2 uðuÞUðu kðbÞÞdu:
ð13Þ
ð14Þ
0
Lemma 6. The function Q(b, b) satisfies the relations Qðb; bÞ ¼ Qðb; bÞ;
ð15Þ
Qðb; bÞ ¼ 1 2UðbÞ þ Qðb; bÞ;
ð16Þ
2
3
Qðb; 0Þ ¼ U ðbÞ;
ð17Þ
Qðb; 1Þ ¼ UðbÞ:
ð18Þ
If b = ±1 the integral in (10) should be understood in the sense of its limit value.
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
1467
Eq. (16) states that the solution b(p, Q) of the calibration problem with input data (p, Q) satisfies the equality bðp; QÞ ¼ bð1 p; Q þ 1 2pÞ:
ð19Þ
Therefore, b(p, Q) is also a solution of the calibration problem with the input data (1 p, 1 2p + Q). The next lemma allows us to find the function Q(b, b) when b = 0. This statement plays the key role in the construction of the approximation of the function Q(b, b), considered in Section 4.1. Lemma 7. For b = 0, the function Q(b, b) satisfies the relation sffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 1 þ b2 p arccosðb2 Þ : ¼ Qð0; bÞ ¼ arctan 2 p 2p 1b
ð20Þ
Let us now prove Theorem 3. Consider Eq. (14). The first term in the right hand side of (14) is given by formula (20). Denote Gðb; cÞ ¼
pffiffiffiffiffiffiffiffi Z 2=p
b
2 =2
eu
c 2 R1 :
Uðu cÞdu;
0
Using this notation, Eq. (14) can be written as Qðb; bÞ ¼ Qð0; bÞ þ Gðb; kðbÞÞ:
ð21Þ
The function G(b, c) satisfies the relation Z b Z Gðb; cÞ ¼ 2 uðuÞUðu cÞdu ¼ 2 0
b
Uðu cÞdUðuÞ:
0
The partial derivative oGðb;cÞ is given by the integral oc Z b oGðb; cÞ 1 1 2 2 ¼2 : uðc uÞuðuÞu du ¼ ð1 eb ð1þc Þ=2 Þ 2 oc p c þ1 0 Since G(b, 0) = U(b) 1/2, we find 1 1 1 Gðb; cÞ ¼ UðbÞ þ arctan c 2 p p
Z
c 0
2
2
eb ðu þ1Þ=2 du: u2 þ 1
ð22Þ
Then from (21), (22) and Lemma 7 we derive4 that 1 1 1 1 arctanð1=kðbÞÞ þ UðbÞ þ arctan kðbÞ p 2 p p Z 1 kðbÞ expðb2 ðu2 þ 1Þ=2Þ ¼p du: p 0 u2 þ 1
Qðb; bÞ ¼
Z
kðbÞ 0
2
2
eb ðu þ1Þ=2 du u2 þ 1
Finally, taking the partial derivative oQðb;bÞ in Eq. (11) we obtain formula (12). The proof of Theorem 3 is now ob complete. h It follows from Eq. (11) that Z 2 1 b z b pffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffi e 1þz2 dz: Qðb; bÞ ¼ U2 ðbÞ þ ð23Þ p 0 1 z2 1 þ z2 Formula (23) can be used for efficient numerical integration of the second moment Q(b, b). In particular, formula (23) is very efficient for small b. In turn, the integral representation in Theorem 3 is very efficient as b 1. Indeed, in this case the upper limit in the integral in (11) is determined by k(b) and k(b) 0. 4
We use the elementary identity arctanw + arctan(1/w) = p/2 in this derivation.
1468
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
3.2. Solvability of the calibration problem Suppose that the empirical mean, p, of the default probabilities and the empirical variance, r2 satisfy the relations 0 < p < 1;
0 6 r2 6 p ð1 pÞ:
ð24Þ
Theorem 8. The calibration problem (6), (7) has a unique nonnegative solution under the condition (24).
Proof. If the empirical variance satisfies (24), the empirical second moment, Q, satisfies the inequalities p2 6 Q 6 p: Since the function Q(b,b) is a continuous, monotonically increasing function of b (see Corollary 4), Q(b, 0) = p2 and Q(b,1) = p, there exists a unique solution of Eq. (7) satisfying the condition 0 6 b 6 1. If Q 62 [p2, p] the calibration problem has no solution. h 4. Approximations of the second moment Q(b, b) 4.1. ‘‘Homotopy approximation’’ of Q(b,b) Now we are ready to construct an approximation of the function Q(b, b). Consider the function QH ðb; bÞ ¼ ð1 kðbÞÞ UðbÞ þ kðbÞ U2 ðbÞ;
ð25Þ
where kðbÞ ¼ p2 arccosðb2 Þ. The function QH(b, b) provides a continuous mapping, QH: b ! QH(b, b), of the unit interval into the subspace of the space of continuous functions such that QH ðb; 0Þ ¼ Qðb; 0Þ ¼ U2 ðbÞ; QH ðb; 1Þ ¼ Qðb; 1Þ ¼ UðbÞ; and for b = 0, QH(0, b) = Q(0, b). The mapping QH is a homotopy in the space of continuous functions. For this reason, it is natural to call QH(b, b) the homotopy approximation of the second moment. From (25) it follows that the approximate solution to the calibration problem satisfies the relation sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p pQ H : ð26Þ b ðp; QÞ ¼ cos 2 p p2 We also note that the function bH is a unique solution to the calibration problem for p = 1/2 (b = 0). Approximation (26) provides an acceptable accuracy for all p 2 [0.3, 0.7]. In this case, the relative error of approximation of Q does not exceed 1%; the relative error of approximate solution of b is smaller than 3%. Even for p = 0.2 the homotopy approximation provides relative error of approximation e = 5% (see Fig. 2, subplot A). For small p this approximation has significant relative error (Fig. 2, subplot C). In this case we construct more complex approximation using asymptotics of Q(b, b) as b ! 0 and b ! 1. 4.2. Approximation for b 0 Eq. (23) implies that if b is close to 0, the function Q(b, b) can be approximated by the integral Z 2 Z b 1 b b2 x eb x 2 ffiffiffiffiffiffiffiffiffiffiffiffi ffi p pffiffiffiffiffiffiffiffiffiffiffiffiffi dx: Qðb; bÞ ¼ p2 þ exp þ dx p 2 p 0 1þx p 0 1 x4 1 x4
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476 A: p= 0.2
B: p= 0.2
1
0.7 0.6 0.5
0.6 β
β(p,Q)
0.8
0.4
0.4 0.3
β(p, Q) H β (p, Q)
0.2 0
1469
0
0.05
0.1 Q
0.15
β(p,Q) β (p,Q) 0 β1(p, Q)
0.2 0.1 0.2
0.04
0.05
0.06
0.07
Q
C: p= 0.033
D: p= 0.033
0.8
0.7
0.6
0.6 β
β(p,Q)
1
0.4
0.5 β(p, Q) βH(p, Q)
0.2 0
0
0.01
0.02 Q
0.03
β(p,Q) β0(p,Q) β1(p, Q)
0.4 0.04
2
4
6 Q
8
10 3
x 10
Fig. 2. Approximations of b(p, Q).
Then we obtain Qðb; bÞ p2 þ
1 b2 e arcsin b2 : 2p
The latter relation implies the following approximation for the solution of the calibration problem: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 bðp; QÞ b0 ðp; QÞ ¼ sinð2p ðQ p2 Þ eb Þ:
ð27Þ
Fig. 2, subplot B, displays the function b0(p, Q). If Q is close to p2 the function b0(p, Q) provides a very good approximation of b(p, Q). 4.3. Approximation for b 1 If b is close to 1, then k(b) is close to 0. Then, from the integral representation (11) we derive Z z expðb2 =2Þ expðb2 t2 =2Þ dt; Qðb; bÞ p p 0
ð28Þ
where z = arctan(k(b)). The integral in (28) is equal Z z pffiffiffiffiffiffi expðb2 t2 =2Þ dt ¼ 2p b1 ðUðbzÞ 1=2Þ: 0
Then, from (28) we derive the asymptotic approximation for the solution of the calibration problem for b 1: bðp; QÞ b1 ðp; QÞ ¼ kðtanð/H =bÞÞ;
ð29Þ
1470
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
where rffiffiffi
p / ¼U jbj exp b2 =2 ðp QÞ þ 0:5 2 H
1
and
rffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 kðxÞ ¼ : 1 þ x2
Fig. 2, subplot D, displays the function b1(p, Q). If Q is close to p the function b1(p ,Q) provides a very good approximation of b(p, Q). 5. Calibration procedure In this Section, we construct an approximate solution, bw(p, Q), of the calibration problem using the approximations bH(p, Q), b0(p, Q) and b1(p, Q). These approximations cover three different areas of values of the parameters p and Q. The function bH provides a very accurate approximation of the function b(p, Q) as p is close5 to 1/2. If Q p2 (b is close to 0), the function b0 approximates the solution of the calibration problem. If Q p (b is close to 1) then the function b1(p, Q) approximates the function b(p, Q). Thus, we have to find a sewing function that provides a smooth transition from one branch of the solution to another. We construct the sewing function using a polynomial regression of the variables sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p pQ H bH ¼ b ðp; QÞ ¼ cos ; 2 p p2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ffi 2 b0 ¼ b0 ðp; QÞ ¼ sin 2p ðQ p2 Þ eb ; b1 ¼ b1 ðp; QÞ ¼ kðtanð/H =bÞÞ; where rffiffiffi
p jbj exp b2 =2 ðp QÞ þ 0:5 / ¼U 2 H
1
and
rffiffiffiffiffiffiffiffiffiffiffiffiffi 1 x2 kðxÞ ¼ : 1 þ x2
The regression coefficients are determined separately in each of the intervals, ½^pðiÞ ; ^pðiþ1Þ , that form a cover of the unit interval, N [ ½^ pðiÞ ; ^ pðiþ1Þ ¼ ½0; 1 i¼1
as a solution to the optimization problem kP2 ðb0 ; b1 ; bH Þ bðp; QÞkL2 ! min;
p0 2 ½^ pðiÞ ; ^pðiþ1Þ ;
where P2 ðx; y; zÞ is a polynomial of degree 2 of the variables x, y and z. Given the empirical default probability, p, and the empirical second moment, Q, we compute b and b as follows: Computation of the parameter b: b = U1(p). Computation of b: If p > 0.5 then p 0 = 1 p and Q 0 = Q + 1 2p otherwise p 0 = p and Q 0 = Q. If 0.3 6 p 0 6 0.5 then bw = bH 5
More precisely, if p 2 [0.3, 0.7] the error of approximation is smaller than 1%.
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
1471
If 0.1 6 p 0 < 0.3, compute bL = a1 Æ b0 + a2 Æ b1 + a3 Æ bH where a1 = 0.3882, a2 = 0.2725 and a3 = 0.3353. if bL > 0:7;
then bH ¼ b1 ;
if 0:7 P bL > 0:3;
then bH ¼ bL ;
else bH ¼ b0 : If 0.05 6 p 0 < 0.1 then if b1 > 0:7;
then bH ¼ b1 ;
elseif b0 < 0:3;
then bH ¼ b0 ;
else bH ¼ bL : If 0.01 6 p 0 < 0.05, compute bQ ¼ a1 b0 þ a2 b1 þ a3 bH þ a4 b0 b1 þ a5 b0 bH þ a6 b1 bH þ a7 b20 þ a8 b21 þ a9 b2H and the coefficients ai are given in the following table a1 = 0.5285 a4 = 1.4945 a7 = 0.0381 if b1 < bQ ;
a2 = 0.36768 a5 = 1.9464 a8 = 0.6292
a3 = 1.7249 a6 = 1.7619 a9 = 1.6296
then bH ¼ b1 ;
elseif bQ > 0:2;
then bH ¼ bQ ;
else bH ¼ b0 : If p 0 < 0.01 we compute the following two quadratic functions bQH ¼ a1 b0 þ a2 b1 þ a3 bH þ a4 b0 b1 þ a5 b0 bH þ a6 b1 bH þ a7 b20 þ a8 b21 þ a9 b2H ; bQL ¼ c1 b1 þ c2 bH þ c3 b1 bH þ c4 b21 þ c5 b2H ; where the coefficients ai are given by a1 = 0.2778 a4 = 1.29100 a7 = 0.3053
a2 = 0.2341 a5 = 0.58969 a8 = 0.42448
a3 = 1.4750 a6 = 1.1624 a9 = 2.4124
c2 = 1.3554 c4 = 1.0509
c3 = 1.0888 c5 = 0.48947
and the coefficients ci are given by c1 = 0.0025767
if bQL P 0:8;
then bH ¼ b1 ;
elseif Q0 P ð0:25 expðb2 Þ þ p2 Þ; elseif bQH P 0:157;
then bH ¼ bQL ;
then bH ¼ bQH ;
else bH ¼ b0 : 6. Numerical results and performance In this Section, we consider the approximate solution to the calibration problem for a small value of the probability of default, p = 0.5%. In Fig. 3 we display the main elements of the approximate solution: the
1472
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476 β~1
β~0
1
0.2
0.15
β(p,Q), β0(p,Q)
β(p, Q), β1(p,Q)
0.8 β1(p, Q)
0.6
β(p,Q)
0.4
0.1
0.05
0.2
β(p, Q) β (p,Q) 0
0
0
1
2 Q
3
4
2.6
2.8
3
3.2
Q
3
x 10
3.4 x 10
5
3.5 3 Relative error (%)
0.85
β
0.8 0.75 True β Approximate β
0.7 0.6
0.8
1
1.2 Q
1.4
2.5 2 1.5 1 0.5 0
1.6
0
0.2
3
x 10
0.4
0.6
0.8
1
β
Fig. 3. Approximation for b(p, Q).
approximating functions b0(p, Q) and b1(p, Q). The function representing the homotopy approximations was already shown in Fig. 2. The first and the second subplots display the approximations b0(p, Q) and b1(p, Q) working for Q p2 and Q p correspondingly. The third subplot represents the resulting approximation obtained by sewing the functions b0(p, Q) and b1(p, Q). The relative error of the approximate solution is shown in the last subplot. The maximal error is attained at the point Qc where the approximate solution switches from the branch corresponding to b0 to the polynomial regression branch linking b0 and b1. In this case, the relative error of the approximation is smaller than 4%. The calibration procedure was tested using exhaustive search over the set of parameters p and b in the area 104 6 p 6 0.5, 104 6 b 6 1 104. We create a grid in the space of parameters on which we study the error of the approximation. First, we computed the value of the second moment Qp(b) = Q(p, b) considered as a function of the parameter b given a value of p. The integral that determines the second moment is computed with the accuracy 1010. After that we apply the calibration procedure to find an approximate value of bw(p, Q). The maximal relative error, e(p), for a given p, p > 0 of the approximation is estimated as eH ¼ maxp eðpÞ, where eðpÞ ¼ max Q>p2
jbðp; QÞ bH ðp; QÞj : bðp; QÞ
Our numerical experiment shows that the relative error of the approximation ew < 0.05. To analyze performance improvement, we measure the computation time of the calibration procedure proposed in Section 5 and the time of the calibration algorithm NINS based on numerical integration and numerical solution of the nonlinear equation ~2 ðp; bÞ r2 ¼ 0; r
ð30Þ
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
1473
Performance ratio
450 400
Perf. ratio
350 300 250 200 150 100 50
0
0.05
0.1
0.15 0.2 0.25 0.3 Default probability
0.35
0.4
0.45
Fig. 4. Performance ratio R = Ts/Ta.
~2 ðp; bÞ ¼ Qðb; bÞ p2 and r is the empirical standard deviation of the default probability. The NINS where r algorithm uses a 30-node Gaussian quadrature (see [10]), Qðb; bÞ
N X i¼1
2
xi U
b bxi ; a
N ¼ 30;
~2 ðp; bÞ. The corresponding summation weights, xi, and the abscissas, xi, are for computation of the variance r calculated by the Gauss–Hermite algorithm. The NINS algorithm finds the positive root of Eq. (30) using the Brent method [1,10]. Then the performance ratio is defined as R = Ts/Ta, where Ta is the time required for the solution of the calibration problem using our approximation and Ts is the computation time of the calibration procedure NINS. We measure the computation time of 100 calibrations corresponding to the values of the parameter b for a given value of the default probability p. The results are shown in Fig. 4. One can see that for the range of values of the default probability p 2 [0.005, 0.5] the performance ratio R > 100. The smallest value is achieved at a very small value of the empirical default probability p = 104. In this case R = 60. Thus, the calibration procedure suggested in this paper demonstrates significant performance advantage and acceptable accuracy.
7. Conclusion We constructed an approximate solution of the calibration problem for the Merton–Vasicek model sewing together three approximations of the second moment of the conditional default probability. This solution allows one to compute the parameters of the model with relative error smaller than 5% even for very small values of the second moment (less than 107). Our solution demonstrates significant performance advantage with respect to the calibration procedure NINS. Our calibration procedure is more than 100 times faster than the calibration algorithm NINS for p > 0.5%.
Acknowledgements The authors thank the referees for the comments and suggestions that helped to improve the paper. We are also very grateful to David Saunders and Ian Iscoe for many useful discussions.
1474
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
Appendix A. Proof of auxiliary results Proof of Proposition 2 Let us compute theRexpectation EUða n þ cÞ, where n Nð0; 1Þ is a standard normal random variable. We 1 have EUða n þ cÞ ¼ 1 Uða x þ cÞdUðxÞ. The integral in this relation can be represented as the probability of the event {f < an + c}, where the random variables f and n are independent and f has a standard normal distribution. Then, from this probabilistic interpretation we obtain that EUða n þ cÞ ¼ Prff < an þ cg ¼ Prff an < cg:
ðA:1Þ pffiffiffiffiffiffiffiffiffiffiffiffiffi The random variable f an has a normal distribution with the mean 0 and the standard deviation 1 þ a2 . Therefore, c Prff an < cg ¼ U pffiffiffiffiffiffiffiffiffiffiffiffiffi : ðA:2Þ 1 þ a2 From (A.1) and (A.2) we obtain (9). Proof of Lemma 5 Let us compute the partial derivative of the function Q(b, b) with respect to the variable b. We have Z 1 oQðb; bÞ o 2 b bx ¼ U dUðxÞ; ob ob a qffiffiffiffiffiffiffiffiffiffiffiffiffi 1 where a ¼ 1 b2 . It is not difficult to see that one can interchange the operators of differentiation and integration in the latter formula. Then we obtain Z oQðb; bÞ 2 1 b bx b bx ¼ U u uðxÞdx ob a 1 a a Z 1 1 b bx 1 b2 b2 þ x2 2bbx þ b2 a2 U exp dx ¼ pa 1 a 2 a2 rffiffiffi Z 2 2 2 b2 =2 1 b bx eðxbbÞ =ð2a Þ pffiffiffiffiffiffi U ¼ e dx: ðA:3Þ p a 2pa 1 Let n Nðbb; aÞ be a random variable having a normal distribution with the expectation bb and the standard deviation a. Denote the distribution function of n by Un(x). Then formula (A.3) can be represented as follows: Z oQðb; bÞ pffiffiffiffiffiffiffiffi b2 =2 1 b bx ¼ 2=p e U ðA:4Þ dUn ðxÞ: ob a 1 The integral in the latter equation satisfies the relation Z 1 b bx U dUn ðxÞ ¼ Pðag þ bn < bÞ; a 1
ðA:5Þ
where g Nð0; 1Þ. The random variable n can be represented in the form n = bb + a Æ f, where f has the standard normal distribution. Then from (A.5) we obtain using Proposition 2 Z 1 qffiffiffiffiffiffiffiffiffiffiffiffiffi b bx U dUn ðxÞ ¼ Pðag þ b ðbb þ afÞ < bÞ ¼ Pðg < ba bfÞ ¼ Uðba= 1 þ b2 Þ ¼ UðbkðbÞÞ; a 1 ðA:6Þ Finally from Eqs. (A.4), (A.5) and (A.6) we find oQðb; bÞ pffiffiffiffiffiffiffiffi b2 =2 ¼ 2=p e Uðb kðbÞÞ: ob
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
1475
Proof of Lemma 6 Let us prove relation (15). We have 0 1 Z 1 B b þ bx C Qðb; bÞ ¼ U2 @qffiffiffiffiffiffiffiffiffiffiffiffiffiA uðxÞdx: 1 1 b2 Changing the integration variable x = t and taking into account that the function u is even, we obtain (15). Formula (16) is proved as follows. For all real x the standard normal cumulative distribution function satisfies the identity U(x) = 1 U(x). Then we have 0 1 0 0 112 Z 1 Z 1 B b bx C B B b þ bx CC Qðb; bÞ ¼ U2 @qffiffiffiffiffiffiffiffiffiffiffiffiffiA uðxÞdx ¼ @1 U@qffiffiffiffiffiffiffiffiffiffiffiffiffiAA uðxÞdx: 1 1 1 b2 1 b2 Then, using Proposition 2 and Eq. (15) we finally obtain (16). Formula (17) follows directly from Eq. (10). To prove (18) notice that qffiffiffiffiffiffiffiffiffiffiffiffiffi 1; if x < b lim Uððb bxÞ= 1 b2 Þ ¼ b!1 0; otherwise: Therefore lim Qðb; bÞ ¼
Z
b!1
b
dUðxÞ ¼ UðbÞ ¼ p:
1
Thus, Lemma 6 is proved. Proof of Lemma 7 qffiffiffiffiffiffiffiffiffiffiffiffiffi Consider the function Q(0, b). Let t ¼ b= 1 b2 . Then the function Q can be written as Z 1 Z 1 Qð0; bÞ WðtÞ ¼ U2 ðtxÞdUðxÞ ¼ U2 ðtxÞdUðxÞ: 1
1
Let us compute the derivative of W(t). We have 2 Z 1 Z dWðtÞ 1 1 x 2 ¼2 UðtxÞ uðtxÞ x uðxÞdx ¼ UðtxÞ exp ðt þ 1Þ x dx; dt p 1 2 1 and integrating by parts we obtain 2 Z 1 dW 1 t x 1 t 1 2 ¼ 2 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi : exp ðt þ 1Þ uðtxÞdx ¼ 2 dt p t þ 1 1 p t þ1 2 2t2 þ 1 Therefore, 1 p
Z
t
s 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ds: 2 þ1 2 þ1 s 2s 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Changing the variable, u ¼ 2s2 þ 1, and taking into account that W(0) = 1/4 we obtain pffiffiffiffiffiffiffiffi Z 2t2 þ1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 1 1 WðtÞ ¼ du ¼ ðarctan 2t2 þ 1 p=4Þ: 2 4 p 1 u þ1 p WðtÞ Wð0Þ ¼
1476
A. Kreinin, A. Nagi / European Journal of Operational Research 185 (2008) 1462–1476
pffiffiffiffiffiffiffiffiffiffiffiffiffiffi Thus, WðtÞ ¼ p1 arctan 2t2 þ 1, and we derive that sffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 1 þ b2 Qð0; bÞ ¼ arctan ; p 1 b2 and after some algebraic transformations we finally obtain (20). Lemma 7 is thus proved. References [1] R.P. Brent, Algorithms for Minimization Without Derivatives, Prentice-Hall, Englewood Cliffs, NJ, 1973. [2] Credit Suisse Financial Products, CreditRisk+—A Credit Risk Management Framework 1997. [3] Klaus Du¨llmann, Harald Scheule, Determinants of the asset correlations of German corporations and implications for regulatory capital, Working paper, University of Regensburg, Department of Statistics, 2003. [4] Michael Gordy, A comparive anatomy of credit risk models, Journal of Banking and Finance 24 (2000) 119–149. [5] Michael Gordy, Erik Heitfield, Estimating default correlations from short panels of credit rating performance data. Working paper, Board of Governors of the Federal Reserve System, 2002. [6] Ian Iscoe, Alexander Kreinin, Default boundary problem, Algorithmics Inc., Research Paper Series, 1999. [7] Ian Iscoe, Alexander Kreinin, Dan Rosen, An integrated market and credit risk portfolio model, Algo Research Quarterly 2 (1) (1999) 21–37. [8] H. Ugur Koyluoglu, Andrew Hickman, A generalized framework for credit risk portfolio models, Risk (2000). [9] Robert Merton, On the pricing of corporate debt: the risk structure of interest rate, Journal of Finance 29 (1974) 449–470. [10] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C, Cambridge University Press, Cambridge, UK, 1993. [11] Rafael Santander, Portfolio Credit Risk Engine, Algorithmics Inc., Technical Report, 2002. [12] Philipp Scho¨nbucher, Credit Derivatives Pricing Models, John Wiley & Sons, 2003. [13] Oldrich Vasicek, Limiting loan loss probability distribution, Technical Document, KMV Corporation, San Francisco, 1991. [14] Oldrich Vasicek, A series expansion for the bivariate normal integral, Technical Document, KMV Corporation, San Francisco, 1998.