Stochastic Local Intensity Loss Models with Interacting Particle Systems

arXiv:1302.2009v1 [math.PR] 8 Feb 2013

Aur´elien Alfonsi ∗ Universit´e Paris-Est, CERMICS, Project Team MathRisk ENPC-INRIA-UMLV, Ecole des Ponts, 6-8 avenue Blaise Pascal, 77455 Marne-la-Vall´ee, France C´eline Labart † Laboratoire de Math´ematiques, CNRS UMR 5127 Universit´e de Savoie, Campus Scientifique 73376 Le Bourget du Lac, France J´erˆome Lelong ‡ Univ. Grenoble Alpes, Laboratoire Jean Kuntzmann, BP 53, 38041 Grenoble, Cedex 09, France February 11, 2013

Abstract It is well-known from the work of Sch¨ onbucher (2005) that the marginal laws of a loss process can be matched by a unit increasing time inhomogeneous Markov process, whose deterministic jump intensity is called local intensity. The Stochastic Local Intensity (SLI) models such as the one proposed by Arnsdorf and Halperin (2008) allow to get a stochastic jump intensity while keeping the same marginal laws. These models involve a non-linear SDE with jumps. The first contribution of this paper is to prove the existence and uniqueness of such processes. This is made by means of an interacting particle system, whose convergence rate towards the non-linear SDE is analyzed. Second, this approach provides a powerful way to compute pathwise expectations with the SLI model: we show that the computational cost is roughly the same as a crude Monte-Carlo algorithm for standard SDEs. Keywords : Stochastic local intensity model, Interacting particle systems, Loss modelling, Credit derivatives, Monte-Carlo Algorithm, Fokker-Planck equation, Martingale problem. AMS classification (2010): 91G40, 91G60, 60H35

Acknowledgments: We would like to thank Benjamin Jourdain (Universit´e Paris-Est) for fruitful discussions having stimulated this research. We would like to thank INRIA Paris Rocquencourt for allowing us to run our numerical tests on the convergence rate on their RIOC cluster. A. Alfonsi acknowledges the support of the “Chaire Risques Financiers” of Fondation du Risque. C. Labart and J. Lelong acknowledge the support of the Finance for Energy Market Research Centre, www.fime-lab.org. ∗

Email:[email protected] Email:[email protected] ‡ Email:[email protected]

1

1

Introduction

In equity modeling, a major concern is to get a model that fits option data. It is well-known from the work of Dupire (1994) that basically European options can be exactly calibrated by using local volatility models σ(t, x). However, local volatility models are known to have some inadequacy to describe real markets. To get richer dynamics, Stochastic Local Volatility (SLV) models have been introduced (see Alexander and Nogueira (2004) or Piterbarg (2006)) and consider the following dynamics for the stock under a risk-neutral probability: dSt = rSt dt + f (Yt )η(t, St )St dWt , where Yt is an adapted stochastic process. Typically, (Yt , t ≥ 0) is assumed to solve an autonomous one dimensional SDE whose Brownian motion may be correlated with W . From the work of Gy¨ongy (1986), we know that under mild assumptions, the following choice η(t, x) = p

σ(t, x) E[f (Yt )2 |St = x]

ensures that St has the same marginal laws as the local volatility model with σ(t, x), which automatically gives the calibration to European option prices. This leads to the following non-linear SDE f (Yt ) σ(t, St )St dWt . dSt = rSt dt + p E[f (Yt )2 |St ]

Here, we stress that the law of (Yt , St ) steps into the diffusion term. Unless for trivial choices of Yt and despite some attempts (Abergel and Tachet (2010)), getting the existence and uniqueness of solutions for this kind of SDE remains an open problem. Also, from a numerical perspective, the simulation of SLV models is not easy, precisely because of the computation of the conditional expectation. In this paper, we propose to tackle a very analogous problem arising in credit risk modeling. In all the paper, we will work under a risk-neutral probabilistic filtered space (Ω, (Ft )t≥0 , F, P). As usual, Ft is the σ-field describing all the events that can occur before time t and F describes all the events. We consider M ∈ N∗ defaultable entities (for example, M = 125 for the iTraxx). We assume that all the recovery rates are deterministic and equal to 1 − LGD for all the firms within the basket. The loss process (Lt , t ≥ 0) is given by Lt =

LGD Xt , ∀t ≥ 0 M

where Xt is the number of defaults up to time t. Clearly, X takes values in LM := {0, . . . , M }. Thanks to the assumption of deterministic recovery rates, and under the assumption of deterministic short interest rates, it is well known that CDO tranche prices only depend on the marginal laws of the loss process (see for example Remark 3.2.1 in Alfonsi (2011)). Let us assume then for a while that we have found marginal laws (P (Xt = k), k ∈ LM )t∈[0,T ] which perfectly fit CDO tranche prices up to maturity T > 0. Then, under some mild assumptions, we know from Sch¨ onbucher (2005) that there exists a non-homogeneous Markov chain with only unit increments which exactly matches these marginal laws. Somehow, this result plays the same role for the loss as Dupire’s result for the stock. The loss model obtained with a non-homogeneous unit-increasing Markov chain is known in the literature as local intensity model. It is fully described by the local intensity λ : 2

R+ × LM → R+ , which gives the instantaneous rate of having one more default in the basket. In the sequel, we will assume that the local intensity λ : R+ × LM → R+ has been calibrated to market data and perfectly matches Index and CDO tranche prices. We make the following assumptions: • ∀x ∈ LM , t ∈ R 7→ λ(t, x) is a c`adl` ag function, • ∀t ≥ 0, λ(t, M ) = 0. In particular, we have λ = maxx∈LM supt∈[0,T ] λ(t, x) < ∞. In this setting, the instantaneous jump rate at time t from x to x + 1 is given by λ(t−, x). Thus, the local intensity model corresponds to a time-inhomogeneous Markov chain making unit jumps with this rate. One may like however get richer dynamics than the ones given by the local intensity model. Then, we can proceed in the same way as for the stochastic volatility models. Let us consider (Yt )t≥0 a general (Ft )-adapted c` adl` ag real process, a function f : R → R+ and a function η : R+ × LM → R+ satisfying the same assumptions as λ. We assume that the default counting process (Xt , t ≥ 0) has jumps of size 1 with the rate η(t−, Xt− )f (Yt− ). By analogy with the equity, we name this kind of model a Stochastic Local Intensity (SLI) model. Then, it is known (see Cont and Minca (2013)) that the local intensity model with the Local Intensity (LI) η(t−, x)E[f (Yt− )|Xt− = x] has the same marginal laws as Xt . Thus, the SLI model will be automatically calibrated to CDO tranche prices if one takes: ∀t > 0, x ∈ LM , η(t−, x) =

λ(t−, x) . E[f (Yt− )|Xt− = x]

This approach has been used in the literature by Arnsdorf and Halperin (2008), and in a slightly different way by Lopatin and Misirpashaev (2008). However, up to our knowledge there is no proof in the literature of the existence nor uniqueness of such a dynamics. The first scope of this paper is to solve this problem. At this stage, we need to make our framework precise. We assume through the paper that: f : R 7→ R is continuous, s.t. ∀x ∈ R, 0 < f ≤ f (x) ≤ f < ∞.

(1.1)

We assume that the probability space (Ω, F, P) contains a standard Brownian motion (Wt , t ≥ 0), a sequence of independent uniform random variables (U k )k∈N , and a sequence (E n )n∈N Pk k n of independent exponential random variables with parameter λf n=1 E f . We set T =

for k ∈ N∗ . The random variables (T k , U k )k∈N∗ will enable us to define a non homogeneous Poisson point process with jump intensity η(t−, Xt− )f (Yt− ). We are interested in studying the following two problems in which we assume that Yt is a process with values either in N (discrete case) or in R (continuous case). In the discrete case, we are interested in finding a predictable process (Xt , Yt )t≥0 such that  P ( )  X = x + k t 0  k,T k ≤t 1  f f (YT k − )λ(T −,XT k − )  U k≤  E[f (Y )|X ] λf

T k−

T k−

 Y0 = y0 , and for each k ≥ 0, (Yt , t ∈ [Tk , Tk+1 )) is a continuous time Markov chain    with transition rate µXt = µXTk . ij

ij

3

(1.2)

For the sake of simplicity, we consider in the discrete setting that X and Y do not jump together almost surely. In the continuous case, the corresponding problem is to solve the following stochastic differential equation:  P ( )  Xt = x0 + k,T k ≤t 1 k f f (YT k − )λ(T k −,XT k − ) U ≤ E[f (Y k )|X k ] λf (1.3) T − T −  Rt Rt Rt  Yt = y0 + 0 b(s, Xs , Ys )ds + 0 σ(s, Xs , Ys )dWs + 0 γ(s−, Xs− , Ys− )dXs ,

for x0 ∈ LM , y0 ∈ R and given real functions b, σ and γ. This framework embeds in particular the dynamics suggested by Arnsdorf and Halperin (2008) and Lopatin and Misirpashaev (2008). Under some rather mild hypotheses on µij , b, σ and γ, which will be specified in the corresponding sections, we will show that the above two equations admit a unique solution. In the discrete case, we are able to show that the corresponding Fokker-Planck equation has a unique solution. This can be achieved by writing the Fokker-Planck equation as an ODE which can be studied directly. In the continuous case, this approach can hardly been extended: the Fokker-Planck equation leads to a non trivial PDE. Instead, we solve this problem by introducing an interacting particle system. This technique is known to be powerful for this type of non linear problems (see Sznitman (1991) or M´el´eard (1996)). The second scope of this paper is to provide a way to compute prices under SLI models. Indeed, interacting particle systems are not only theoretical tools to prove existence and uniqueness results for such equations. They give a very smart way to simulate these processes, therefore enabling us to run Monte-Carlo algorithms. This approach has been recently used by Guyon and Henry-Labord`ere (2011) for Stochastic Local Volatility models. For the Stochastic Local Intensity models considered in this work, the conditional expectation is much simpler to handle. This enables us to get theoretical results on the convergence and also simplifies the implementation. In fact, we show in our case under some assumptions that the rate of convergence to estimate expectations is in O(1/N α ) for any α < 1/2, where N is the number of particles. On our numerical experiments, we even observe on several examples a convergence which is similar to the one of the Central Limit Theorem, which is rather usual for Interacting Particle Systems. Besides, we show that we can simulate the interacting particle system with a computational cost in O(DN ), where D is the number of time steps for the discretization of the SDE on Y . Thus, the computational cost is roughly the same as a crude Monte-Carlo algorithm for standard SDEs with N samples. The paper is organized as follows. First, we study the case where Y has discrete values; this framework enables us to settle the problem and solve it by rather elementary tools. This part is independent from the rest of the paper. Second, by means of a particle system approach, we investigate the case where Y is real valued jump diffusion. Finally, we carry out numerical simulations highlighting the relevance of the particle system technique to compute pathwise expectation of the Process (1.3).

2

The SLI model when Y takes discrete values

The goal of this section is to prove the existence of a process (Xt , Yt )t≥0 satisfying (1.2). Unlike the continuous case (1.3), we can get this result by elementary means, without resorting to 4

an interacting particle system. To do so, we write the Fokker-Planck equation associated to the process (Xt , Yt )t≥0 , which should be satisfied by P(Xt = i, Yt = j) for (i, j) ∈ LM × N. We have  P λ(t,i−1) i  ∂t p(t, i, j) = f (j)p(t, i − 1, j)  k6=j µkj p(t, i, k) + 1{i≥1} ϕp (t,i−1)      − ϕλ(t,i) f (j)1{i≤M −1} − µijj p(t, i, j) (E) p (t,i)   p(0, i, j) = 0 ∀(i, j) 6= (x0 , y0 ),    p(0, x , y ) = 1. 0 0 where p is a function from R+ × LM × N to R and P∞ j=0 f (j)p(t, i, j) . ϕp (t, i) = P∞ j=0 p(t, i, j)

If we manage to prove that the Fokker-Planck equation admits a unique solution p such that ∀i, j ∈ LM × N, ∀t ≥ 0, p(t, i, j) ≥ 0,

(2.1)

∀t ≥ 0

(2.2)

M X ∞ X

p(t, i, j) = 1,

i=0 j=0

then we will get that the law of a process (Xt , Yt )t≥0 satisfying (1.2) is unique. Besides, we will also get the existence of such a process. It is easy to check that a continuous Markov chain (Xt , Yt )t≥0 starting from (x0 , y0 ) with transition rate matrix µ ˜(i1 ,j1 ),(i2 ,j2) (t) = 1j1 =j2 (1i2 =i1 +1 − 1i2 =i1 )

f (j1 )λ(t, i1 ) + 1i1 =i2 µij11 j2 , i1 , i2 ∈ LM , j1 , j2 ∈ N, ϕp (t, i1 )

where p is the solution of (E) satisfies (1.2). In fact, the Fokker-Planck equation of this process  P λ(t,i−1) i  f (j)q(t, i − 1, j)  ∂t q(t, i, j) = k6=j µkj q(t, i, k) + 1{i≥1} ϕp (t,i−1)      − ϕλ(t,i) f (j)1{i≤M −1} − µijj q(t, i, j) p (t,i)   = 0 ∀(i, j) 6= (x0 , y0 ),  q(0, i, j)   q(0, x , y ) = 1. 0 0 is linear and clearly solved by p, which gives q ≡ p.

2.1

Assumptions and notations

In this part, we assume that the transition rates satisfy the following hypothesis. Hypothesis 1 The intensity matrices (µkij )i,j≥0 for k ∈ LM satisfy the following conditions: • ∀k ∈ LM , ∀i, j ∈ N × N such that i 6= j µkij ≥ 0, • ∀k ∈ LM , ∀i ∈ N µkii ≤ 0, P P∞ k k k • ∀k ∈ LM , ∀i ∈ N ∞ j=0 µij = 0 (then j=0,j6=i µij = −µii ). 5

Moreover, we assume ∀k ∈ LM , supi∈N |µkii | < ∞. We also introduce specific notations used in the discrete case. • We define the set E as the set of real sequences indexed by i ∈ LM and j ∈ N: E := {u = (uij )0≤i≤M,0≤j : uij ∈ R},

∗ := {u = (ui ) i ∗ E+ := {u = (uij )0≤i≤M,0≤j : uij ∈ R+ } and E+ j 0≤i≤M,0≤j : uj ∈ R+ }. P P∞ i • For u ∈ E, we set |u| = M i=0 j=0 |uj |.

2.2

Solving the Fokker-Planck equation

To be more concise, we rewrite the Fokker-Planck Equation (E) and Conditions (2.1) − (2.2) by using a sequence of functions. Let P := (Pji )0≤i≤M,j≥0 denote a sequence such that each Pji is a function from R+ to R. Solving (E) under Constraints (2.1) − (2.2) boils down to solving ( P ′ (t) = Ψ(t, P (t)), ′ (E ) P (0) = P0 under the constraints ∀t ≥ 0, P (t) ≥ 0 and |P (t)| = 1. The sequence P0 is such that (P0 )ij = 0 for (i, j) 6= (x0 , y0 ) and (P0 )xy00 = 1. Ψ is an application from R+ × E+ → E given by   X λ(t, i − 1) λ(t, i) (Ψ(t, x))ij = µikj xik + 1{i≥1} f (j)xi−1 − 1 f (j) xij {i≤M −1} j ϕ(x, i − 1) ϕ(x, i) where ϕ(x, i) :=

k≥0 P∞ i l=0 f (l)xl P ∞ i . x l=0 l

∗ . When x ∈ E , difficulties may arise Remark 1. Ψ is defined without ambiguity on E+ + zji when for some fixed i, xij = 0 ∀j. In this case, we still have lim = 0 since ∗ ,z→x ϕ(z, i) z∈E+ ∗ . Thus, we can extend Ψ by continuity on E . f ≤ ϕ(z, i) ≤ f for z ∈ E+ +

We aim at solving (E) in the set of summable sequences (compatible with Condition (2.2)) and get the following result.

Theorem 2. Equation (E ′ ) admits a unique solution on R+ satisfying ∀t ≥ 0, P (t) ≥ 0 and |P (t)| = 1. To do so, we first focus on the following differential equations: ( P ′ (t) = Ψ(t, (P (t))+ ), ′′ (E ) P (0) = P0 .

Proposition 3. Equation (E ′′ ) admits a unique solution on R+ . Moreover, the solution satisfies ∀t ≥ 0, P (t) ≥ 0 and |P (t)| = 1.

The proof of this proposition consists in first showing the Lipschitz property which gives the existence and uniqueness of P and then proving that ∀t ≥ 0, P (t) ≥ 0 and |P (t)| = 1. This proof is postponed to Appendix A. Then, the proof of Theorem 2 becomes obvious. The unique solution P of (E ′′ ) clearly solves (E ′ ), and any solution Q of (E ′ ) such that ∀t ≥ 0, Q(t) ≥ 0 and |Q(t)| = 1 also solves (E ′′ ) and thus coincides with P . 6

3

The SLI model when Y is real valued.

3.1

Setting and main results

We are interested in proving the existence of a process (Xt , Yt )t solving the stochastic differential equation (1.3). More precisely, we will consider the following stochastic differential equation,  P ( )  Xt = X0 + k,T k ≤t 1 k f f (YT k − )λ(T k −,XT k − ) U ≤ E[f (Y k )|X k ] λf (3.1) T − T −  Rt Rt Rt  Yt = Y0 + 0 b(s, Xs , Ys )ds + 0 σ(s, Xs , Ys )dWs + 0 γ(s−, Xs− , Ys− )dXs ,

with (possibly) random initial condition (X0 , Y0 ) such that E[|Y0 |m ] < ∞ for any m ∈ N. We will denote in the sequel Linit the probability law of (X0 , Y0 ) under P. To get existence and uniqueness results for (3.1), we will make the following assumption on the coefficients. Hypothesis 2 1. The functions b, σ, γ : R+ × LM × R → R are measurable, with sublinear growth with respect to y: ∀T > 0, ∃CT > 0, ∀t ∈ [0, T ], x ∈ LM , y ∈ R, |b(t, x, y)|+|σ(t, x, y)|+|γ(t, x, y)| ≤ CT (1+|y|). 2. The functions b(t, x, y) and σ(t, x, y) are such that for any x ∈ LM , y0 ∈ R, there exists a unique strong solution for the SDE Z t Z t σ(s, x, Ys )dWs , t ≥ 0. b(s, x, Ys )ds + Yt = y 0 + 0

0

This property holds if we assume for example that: ( |b(t, x, y) − b(t, x, y ′ )| ≤ CT |y − y ′ | p ∀T > 0, ∃CT > 0, ∀t ∈ [0, T ], x ∈ LM , |σ(t, x, y) − σ(t, x, y ′ )| ≤ CT |y − y ′ |.

3. For any x ∈ LM , (t, y) 7→ γ(t, x, y) is c` adl` ag with respect to t and continuous with respect to y, i.e. γ(t, x, y) = lims>t,s→t,z→y γ(s, x, z) and lims 0 and ϕQ (t, x) = f otherwise. Q(Xt = x)

Let Hypothesis 2 hold. Then, for any (x0 , y0 ) ∈ LM × R, there exists a unique strong solution (Xt , Yt )t∈[0,T ] to the following SDE with jumps:  P ) (  Xt = x0 + k,T k ≤t 1 k f f (YT k − )λ(T k −,XT k − ) U ≤ λf ϕQ (T k −,X k ) (3.3) T −  Rt Rt Rt  Yt = y0 + 0 b(s, Xs , Ys )ds + 0 σ(s, Xs , Ys )dWs + 0 γ(s−, Xs− , Ys− )dXs . Proof of Theorem 4. The proof of Theorem 4 is split in three main steps.

• First, we show the existence of a probability measure Q ∈ P(E) solving the Martingale Problem (3.2). This result is obtained by considering the associated interacting particle system: we show that each particle converges in law, and that any probability measure in the support of the limiting law solves the Martingale Problem. This is done in Section 3.3. • Second, we show the uniqueness of the probability measure Q ∈ P(E) solving (3.2). To do so, we introduce a function Ψ : P(E) → P(E) defined as follows. Let (X0 , Y0 ) be a random variable distributed according to Linit under P. Then, we know from Proposition 5 that there exists a unique process (XtQ , YtQ )t∈[0,T ] solving  P Q     Xt = X0 + k,T k ≤t 1 k f f (YTQk − )λ(T k −,XTQk − )  U ≤

λf

Q ϕQ (T k −,X

)

T k−   Y Q = Y + R t b(s, X Q , Y Q )ds + R t σ(s, X Q , Y Q )dW + R t γ(s−, X Q , Y Q )dX Q , s 0 s s s s s s− s− t 0 0 0 (3.4)



8



and we define Ψ(Q) = law((XtQ , YtQ )t∈[0,T ] ) ∈ P(E). As we will see in Section 3.4, Ψ (or more precisely Ψ iterated k-times) is a contraction mapping for the variation norm. Combining this result with the following Lemma gives the uniqueness of the probability measure solving (3.2). Lemma 6. Let Q ∈ P(E). We have Ψ(Q) = Q if and only if Q solves the Martingale Problem (3.2). In this case, (XtQ , YtQ )t∈[0,T ] solves Equation (3.1). • Then, the existence and uniqueness of Q satisfying Ψ(Q) = Q and Lemma 6 automatically give the strong existence and uniqueness of (X, Y ) satisfying (3.1) since we necessarily have E[f (YT k − )|XT k − ] = EQ [f (YT k − )|XT k − ] and thus (Xt , Yt )t∈[0,T ] = (XtQ , YtQ )t∈[0,T ] .

 Proof of Lemma 6. The direct implication is clear since Ψ(Q) = Q gives that E[f (YtQ )|XtQ ] = ϕQ (t, XtQ ). Thus, (XtQ , YtQ )t∈[0,T ] solves (3.1) and in particular Q solves the Martingale Problem (3.2). Conversely, let us assume that Q solves the Martingale Problem (3.2). We know from Proposition 5 that strong uniqueness holds for the SDE (XtQ , YtQ )t∈[0,T ] . From Lepeltier and Marchal (1976, Theorem II13 ), we know that strong uniqueness implies weak uniqueness, which precisely gives Ψ(Q) = Q since Q solves the Martingale Problem (3.2) and therefore the Martingale problem associated with (3.4). In particular, we have E[f (YtQ )|XtQ ] = ϕQ (t, XtQ ), and (XtQ , YtQ )t∈[0,T ] solves (3.1). 

3.2

The interacting particle system

We assume that the probability space (Ω, F, P) carries all the random variables used below. Now, we set up the particle system related to the Martingale Problem (3.2). In the following, N will denote the number of particles, (X0i , Y0i ), i ∈ N are independent random variables following the law Linit under P, and (Wti , t ≥ 0), i ∈ N are independent standard Brownian motions. We build an interacting particle system ((Xti,N , Yti,N ), t ≥ 0)1≤i≤N with the following features. For i = 1, . . . , N , (Xti,N , t ≥ 0) is a Poisson process with intensity: i,N i,N PN λ(t−, Xt− )f (Yt− ) j=1 1{X j,N =X i,N } t− t− , (3.5) PN j,N f (Y )1 j,N i,N t− j=1 {X =X } t−

t−

and Yti,N solves the following equation: Yti,N = Y0i +

Z

t 0

b(s, Xsi,N , Ysi,N )ds +

Z

t 0

σ(s, Xsi,N , Ysi,N )dWsi +

Z

t 0

i,N i,N )dXsi,N . γ(s−, Xs− , Ys−

(3.6)

In fact, we can give an explicit construction of this particle system, which will be useful later and we explain now. Let us consider (U i,k )i,k∈N∗ a sequence of independent uniform 9

variables on [0, 1] and (E i,k )i,k∈N∗ a sequence of independent exponential random variables with parameter λf f . These variables are independent, and independent of the previously defined Brownian motions. We define the times T i,k =

k X

E i,l ,

l=1

and we can order (T i,k , i = 1, . . . , N, k ≥ 1), such that 0 < T i1 ,k1 < · · · < T il ,kl < . . . almost surely. Up to the first jump of X i,N , Yti,N is defined as the unique strong solution of Z t Z t σ(s, X0i , Ysi,N )dWsi . b(s, X0i , Ysi,N )ds + Yti,N = Y0i + 0

0

At time τ = T il ,kl , the process X il ,N makes a jump of size 1 if P l ,N ) N λ(τ −, Xτil−,N )f (Yτi− j=1 1{X j,N =X il ,N } f τ− τ− , U il ,kl ≤ PN j,N λf )1 f (Y il ,N j,N τ− j=1 {X =X } τ−

τ−

l ,N l ,N ) + γ(τ −, Xτil−,N , Yτi− and does not jump otherwise. If a jump occurs, we set Yτil ,N = Yτi− il ,N i ,N l and, up to the next jump of X , we define Yt as the unique strong solution of Z t Z t σ(s, Xτil ,N , Ysil ,N )dWsi , t ≥ τ. b(s, Xτil ,N , Ysil ,N )ds + Ytil ,N = Yτil ,N +

τ

τ

3.3

Existence of a solution to (3.2)

We followPthe analysis carried out by M´el´eard (1996), pages 69 and 70. We denote by µN = N1 N the empirical measure given by the particle system. It is a i=1 δ(X i,N ,Y i,N ) t

t

t∈[0,T ]

random variable taking values in P(E). We denote by π N ∈ P(P(E)) the probability law of µN . For π ∈ P(P(E)), we denote by Z µπ(dµ) ∈ P(E), I(π) = P(E)

the mean of π. Let F : E → R be a bounded function which is continuous with respect to the Skorokhod topology. It induces an application P(E) → R — still denoted by F with an abuse of notation — such that Z F (z)µ(dz). F (µ) = E

Since π N is by definition the probabilityR law of µN , we have by using Fubini’s TheoR N rem E(F (µ )) = P(E) F (µ)π N (dµ) = E F (z)I(π N )(dz). On the other hand, we have P i,N i,N )t∈[0,T ] ). By symmetry, (Xti,N , Yti,N )t∈[0,T ] has the same law F (µN ) = N1 N i=1 F ((Xt , Yt as (Xt1,N , Yt1,N )t∈[0,T ] and we get that: Z 1,N 1,N F (z)I(π N )(dz). (3.7) E[F ((Xt , Yt )t∈[0,T ] )] = E

10

Lemma 7. The sequence (π N )N is tight. The proof of Lemma 7 is postponed to Appendix D. Now, wWe can consider a subsequence π Nk which converges weakly to π ∞ ∈ P(P(E)). Let q ∈ N∗ , 0 ≤ s1 ≤ · · · ≤ sq ≤ s ≤ t, φ ∈ Cb0,2 (R2 , R) and g1 , . . . , gq ∈ Cb (R2 , R) be bounded functions with bounded derivatives for φ. We set for Q ∈ P(E), Z tn λ(u, Xu )f (Yu ) (φ(Xu + 1, Yu + γ(u, Xu , Yu )) − φ(Xu , Yu )) EQ [f (Yu )|Xu ] s q o Y i 1 2 2 gl (Xsq , Ysq ) (3.8) + b(u, Xu , Yu )∂y φ(Xu , Yu ) + σ (u, Xu , Yu )∂y φ(Xu , Yu ) du 2

F (Q) =EQ

h

φ(Xt , Yt ) − φ(Xs , Ys ) −

l=1

We have to check that Q 7→ F (Q) is continuous with respect to Q for the weak convergence. Since f is continuous by Assumption (1.1), we first notice that when Qn converges weakly to Q, EQn [f (Yt )|Xt = x] converges to EQ [f (Yt )|Xt = x] when Q(Xt = x) > 0, unless for an at most countable set of times t depending on Q. Then, following M´el´eard (1996), it comes out that if s, t, s1 , . . . , sq are taken outside a countable set depending on π ∞ , F is π ∞ -a.s. continuous. In this case, we have: Z F (Q)2 π ∞ (dQ). E[F (µNk )2 ] → k→+∞

P(E)

By definition of µN , we have: q N i Y 1 X h i,N i,N ) , where , Y gl (Xsi,N (Mt − Msi,N ) sq q N i=1 l=1 Z tn 1 b(u, Xui,N , Yui,N )∂y φ(Xui,N , Yui,N ) + σ 2 (u, Xui,N , Yui,N )∂y2 φ(Xui,N , Yui,N ) Mti,N = φ(Xti,N , Yti,N ) − 2 0 i,N i,N PN o λ(u, Xu )f (Yu ) j=1 1{X j,N =X i,N } u u i,N i,N i,N i,N i,N i,N [φ(X + 1, Y + γ(u, X , Y )) − φ(X , Y )] du. + PN u u u u u u j,N j=1 f (Yu )1{X j,N =X i,N }

F (µN ) =

u

u

We observe now that [M i,N , M j,N ]t = 0 for i 6= j since, by construction these martingales do not jump together and hW i , W j it = 0. Therefore, we get that: q

i Y 1 h 1,N 2 E[F (µ ) ] = E (Mt1,N − Ms1,N )2 ≤ C/N, gl (Xs1,N , Y ) sq q N N 2

l=1

thanks to the boundedness assumption made on functions gl and φ. It comes out that F (Q) = 0, π ∞ (dQ) almost surely. This holds in fact for any function F given by (3.8), provided that s, t, s1 , . . . , sq are taken outside a countable set depending on π ∞ . Since the process (Xt , Yt ) is c`adl` ag, this is sufficient to show that any measure in the support of π ∞ solves the Martingale Problem (3.2). In particular, we get the following result. Proposition 8. There exists a measure Q ∈ P(E) solving the Martingale Problem (3.2).

11

3.4

Uniqueness of a solution to (3.2)

Let Q ∈ P(E) denote a probability measure solving the Martingale Problem (3.2). We know that such a probability exists thanks to Proposition 8. We want to show that it is indeed unique. To do so, we consider another probability Q ∈ P(E) and study the total variation distance VT (Ψ(Q) − Ψ(Q)) between Q and Q over E. For t ∈ [0, T ], Q, Q ∈ P(E), we denote by Q [0,t] and Q [0,t] their restriction to the paths on the time interval [0, t]. We also set and Q . Vt (Q − Q) the total variation distance between Q [0,t]

[0,t]

Lemma 9. Let τ := inf{t ≥ 0, XtQ 6= XtQ } denote the first time when X Q and X Q do not jump together. We have VT (Ψ(Q) − Ψ(Q)) ≤ 2P(τ ≤ T ). Proof. Let us recall that for any signed measure η on E, the total variation of η is given by V (η) = η + (E) + η − (E), where η = η + − η − is the Hahn-Jordan decomposition of η. Besides, we clearly have 1 V (η) ≤ V˜ (η) ≤ V (η), 2 where V˜ (η) = sup{|η(A)|, A ⊂ E measurable}. We have for any measurable set A of E, P((XtQ , YtQ )t∈[0,T ] 6= (XtQ , YtQ )t∈[0,T ] )

≥ P((XtQ , YtQ )t∈[0,T ] ∈ A, (XtQ , YtQ )t∈[0,T ] 6∈ A)

= P((XtQ , YtQ )t∈[0,T ] ∈ A) − P((XtQ , YtQ )t∈[0,T ] ∈ A, (XtQ , YtQ )t∈[0,T ] ∈ A)

≥ P((XtQ , YtQ )t∈[0,T ] ∈ A) − P((XtQ , YtQ )t∈[0,T ] ∈ A). By taking the supremum over A, we get

P((XtQ , YtQ )t∈[0,T ] 6= (XtQ , YtQ )t∈[0,T ] ) ≥ V˜T (Ψ(Q) − Ψ(Q)), which gives the claim.



Up to time τ , we have (XtQ , YtQ ) = (XtQ , YtQ ), and we get that Z t 1 1 Q Q 1{τ ≤t} − 1{τ >s} λ(s, Xs )f (Ys ) − ds Q Q 0 ϕQ (s, Xs ) ϕQ (s, Xs )

is a martingale. In particular, we have # " 1 dP(τ ≤ t) 1 Q Q − = E 1{τ >t} λ(t, Xt )f (Yt ) ϕQ (t, X Q ) ϕQ (t, X Q ) dt t t i λf h Q E |ϕ (t, XtQ ) − ϕQ (t, XtQ )| . ≤ f2 12

Now,

let us observe that ϕQ (t, x) − ϕQ (t, x)

=

EQ [f (Yt )1{Xt =x} ]−EQ [f (Yt )1{Xt =x} ] Q(Xt =x) Q |ϕ (t, x) −

EQ [f (Yt )|Xt =x] [Q(Xt = x) − Q(Xt = x)]. Thus, we get that Q(Xt =x)  1 |EQ [f (Yt )1{Xt =x} ] − EQ [f (Yt )1{Xt =x} ]| + f |Q(Xt = x) − Q(Xt Q(X =x) t

fore (we use here that Q is invariant, and is the law of (XtQ , YtQ )t≥0 ) X

E[|ϕQ (t, XtQ ) − ϕQ (t, XtQ )|] ≤

x∈LM

+

ϕQ (t, x)| ≤  = x)| , and there-

|EQ [f (Yt )1{Xt =x} ] − EQ [f (Yt )1{Xt =x} ]| + f |Q(Xt = x) − Q(Xt = x)|

≤ 2f Vt (Ψ(Q) − Ψ(Q)). To Pncheck the last inequality, one has to observe that for a simple nonnegative function f (y) = fi 1{Ai } , where Ai ⊂ R are Borel sets, we have |EQ [f (Yt )1{Xt =x} ] − EQ [f (Yt )1{Xt =x} ]| ≤ Pni=1 i=1 fi |Q(Xt = x, Yt ∈ Ai ) − Q(Xt = x, Yt ∈ Ai )| ≤ f Vt (Q − Q). By passing to the limit, this property holds for any bounded measurable (hence continuous) function f . To sum up, we have for t ∈ [0, T ], λf VT (Ψ(Q) − Ψ(Q)) ≤ 2P(τ ≤ T ) ≤ 4 2 f

2

Z

2

T

0

Vs (Q − Q)ds ≤ 4

λf T VT (Q − Q), f2

since V0 (Ψ(Q) − Ψ(Q)) = 0 and s 7→ Vs (Q − Q) is nondecreasing. Let us recall that the set of bounded countably additive measures on E endowed with the total variation norm is a 2

Banach space. If 4 λff2 T < 1, we get by the Banach fixed point theorem that Ψ has a unique fixed point that is necessarily Q. Otherwise, we get by iterating that:

VT (Ψ

(k)

(Q) − Ψ

(k)

λf (Q)) ≤ 4 2 f 



f f2

2

2

T

Z

0

T

Vs (Ψ(k−1) (Q) − Ψ(k−1) (Q))ds ≤

k  2 λ f 4 f2 T k!

VT (Q − Q).

k

< 1 and Ψ(k) is a contraction mapping. Thus, Ψ(k) has a When k is large enough, k! unique fixed point that is necessarily Q. This concludes the proof of Theorem 4. Besides, using the notations of Section 3.3, we get that any convergent subsequence of π N should converge to δQ (dQ). This gives the weak convergence of π N towards δQ (dQ). By (3.7), we get that any particle converges in law towards Q.

3.5

Convergence speed towards Q

Now that we have proved that each particle converges to the invariant probability measure, we are interested in characterizing the speed of convergence of the interacting particle system towards this measure. This question is of practical importance, since one would like to use the following approximation EQ [F ((Xt , Yt )t∈[0,T ] )] ≈

N 1 X F ((Xti,N , Yti,N )t∈[0,T ] ), N i=1

13

(3.9)

and have an estimate of the error involved. First, we need to introduce some additional notations. We consider the same particle system (Xti,N , Yti,N )t≥0 as in Section 3.2, constructed with the random variables T i,k , U i,k and ¯ i , Y¯ i )t≥0 as the unique (Wti , t ≥ 0). With these variables, we construct now the processes (X t t solution of  P   ¯ ti = X ¯i +  X  ¯i ¯i 0 k,T i,k ≤t 1 f (Y )λ(T i,k −,X )  i,k i,k f T − T − U i,k ≤ λf ¯i (3.10)   ϕQ (T i,k −,X ) i,k − T  Rt Rt  Y¯ i = Y¯ i + R t b(s, X i i i i i i i i ¯ . ¯ , Y¯ )dX ¯ , Y¯ )dW + γ(s−, X ¯ , Y¯ )ds + σ(s, X t

0

s

s

0

s

s

0

s

0

s−

s−

s

¯ i , Y¯ i )t∈[0,T ] is the invariant probability law Q, since Ψ(Q) = Q. By construction, the law of (X t t By using the same argument as in Lemma 9, we have: ¯ 1 , Y¯ 1 )t∈[0,T ] ) = 2P(τ 1 ≤ T ), VT (L((Xt1,N , Yt1,N )t∈[0,T ] ) − Q) ≤ 2P((Xt1,N , Yt1,N )t∈[0,T ] 6= (X t t ¯t1 6= X 1,N }. We also set for i = 2, . . . , N , τ i = inf{t ≥ 0, X ¯ti 6= X i,N }. where τ 1 = inf{t ≥ 0, X t t Proposition 10. Let us assume that Linit is such that P(X0 = x) > 0 for any x ∈ LM . Then, there is a constant K > 0 such that K P(τ 1 ≤ T ) ≤ √ . N ¯ 1 and X 1,N may become different at the times T 1,k if Proof. By construction, the processes X U 1,k



1,k ¯1 ¯1 f λ(T −,XT 1,k − )f (YT 1,k − ) 1 ¯ λf ϕQ (T 1,k −,X )

and

U 1,k

T 1,k −

Thus, 1{τ 1 ≤t} −

Rt 0

f

1{τ 1 >s} λf

have

1

P(τ ≤ T ) =



λ(s,X¯ s1 )f (Y¯s1 ) ϕQ (s,X¯ 1 ) − s

Z T  E λf   0 f

>

1,N λ(T 1,k −,X 1,N 1,k )f (Y 1,k )

f λf

T − T − PN j,N f (Y )1 j,N 1,N j=1 T 1,k − {X =X } T 1,k − T 1,k − PN 1 1,N j=1 {X j,N =X } T 1,k − T 1,k −

, or conversely.

λ(s,Xs1,N )f (Ys1,N ) ds is a martingale and we PN j,N f (Ys )1 j,N 1,N j=1 {Xs =Xs } PN 1 j=1 {X j,N =X 1,N } s

¯t1 )f (Y¯t1 ) λ(t, X − 1{τ 1 >t} ¯ 1) ϕQ (t, X t

s

  1,N 1,N λ(t, Xt )f (Yt )   PN j,N dt . )1 j,N 1,N j=1 f (Yt  {Xt =Xt } PN j=1 1{X j,N =X 1,N } t

t

¯t1 , Y¯t1 ) = (X 1,N , Y 1,N ). Therefore, Let us observe that on {τ 1 > t}, we have (X t t   PN Z T 1 + 1 j,N 1 ¯ j=2 {Xt =Xt } 1 1 dt . (3.11) 1{τ 1 >t} − P(τ ≤ T ) ≤ f E  P ¯ 1 ) f (Y 1,N ) + N f (Y j,N )1 j,N ¯ 1 0 ϕQ (t, X t t t j=2 =X } {X t

Now, we study Q 1 − ϕ (t,x)

P 1+ N j=2 1{X j,N =x} t PN for x ∈ LM , and set j,N 1,N f (Yt )+ j=2 f (Yt )1 j,N {X =x} t

q¯t (x) = Q(Xt = x). 14

t

When q¯t (x) > 0, we have   PN q¯ (x) − 1 1 + PN 1 j,N 1 + 1 j,N t j=2 {Xt =x} j=2 {Xt =x} N 1 P ϕQ (t, x) − ≤ j,N [f (Y )1 ] E t f (Yt1,N ) + N f (Y )1 {X =x} j,N t Q t j=2 {Xt =x}   N X 1 1  j,N 1,N  + EQ [f (Yt )1{Xt =x} ] − f (Yt )1{X j,N =x} f (Yt ) + t f EQ [f (Yt )1{Xt =x} ] N j=2   N X 1 1  ≤ 1{X j,N =x}  q¯t (x) − 1+ t f q¯t (x) N j=2   N X 1 1 j,N 1,N f (Yt ) + f (Yt )1{X j,N =x}  . + 2 EQ [f (Yt )1{Xt =x} ] − t N f q¯t (x) j=2

¯ N +1 , Y¯ N +1 )t∈[0,T ] another We analyze these two terms in a similar manner. We introduce (X t t ¯ 1 , Y¯ 1 )t∈[0,T ] , which is independent from all other existing processes. We have: copy of (X t t   N N N +1 X X X 1 1 1 1 q¯t (x) − 1 + 1{X¯ j =x} − 1{X j,N =x} + 1{X j,N =x}  ≤ q¯t (x) − 1{X¯ j =x} + t t t t N N N N j=2 j=2 j=2 N +1 N 1 1 X 1 X ≤ q¯t (x) − (3.12) 1{X¯ j =x} + 1{τ j ≤t} + , t N N N j=2

j=2

¯ j and X j,N may be different only on {τ j ≤ t}. Similarly, we have since X t t

We introduce

  N X j,N E [f (Yt )1{X =x} ] − 1 f (Y 1,N ) +  f (Yt )1{X j,N =x} t t Q t N j=2 N +1 N f X f 1 X ¯j f (Yt )1{X¯ j =x} + 1{τ j ≤t} + . ≤ EQ [f (Yt )1{Xt =x} ] − t N N N j=2 j=2

(3.13)

  N +1 N +1 X X 1 1  1 1 j ¯  . [f (Y )1 ] − A(x) = 1 f ( Y )1 + q ¯ (x) − E j j t t ¯ ¯ {X =x} t t {Xt =x} {Xt =x} q¯t (x) N f Q N j=2 j=2

¯ 1 , Y¯ 1 )t≥0 . From (3.11), (3.12) This is a random function which is independent from (X t t j 1 and (3.13), we obtain by observing that τ and τ have the same law under P:       Z T 1 1 f 1 1 1 ¯ E E[A(Xt )] + 1 + P(τ ≤ T ) ≤ (3.14) ¯ 1 ) + E q¯t (X ¯ 1 ) 1τ 2 ≤t dt. f N q¯t (X 0 t t P q¯t (x) 1 First, we observe that E[ q¯ (X ¯ 1) ] = x, s.t. q¯t (x)>0 q¯t (x) ≤ M + 1. Thanks to the independence t t ¯ 1 = x] = E[A(x)]. On the one hand, we have: ¯ 1 )|X ¯ 1 , E[A(X of A(x) and X t

t

t

15

h E q¯t (x) −

1 N

P N +1 j=2

 i 1{X¯ j =x} ≤ t

s   E q¯t (x) −

On the other hand, we have similarly that:

2 

1 PN +1 ¯ j =x} ) N ( j=2 1{X t

=



q¯t (x)(1−¯ qt (x)) √ . N

    N +1 X 1 j ¯  f (Yt )1{X¯ j =x}   E  EQ [f (Yt )1{Xt =x} ] − t N j=2 v  s 2  r u u 2 N +1 X u 1 1 f ≤ tE EQ [f (Yt )1{Xt =x} ] − f (Y¯tj )1{X¯ j =x} )  ≤ EQ (1{Xt =x} f 2 (Yt )) ≤ q¯t (x). t N N N j=2

Finally, we obtain that:

E[A(x)] ≤



f 1+ f



p

1 √ . q¯t (x) N

By using the tower property of the conditional expectation, we get: √  f M +1 1 ¯ √ , E[A(Xt )] ≤ 1 + f N p √ P since we have x∈LM q¯t (x) ≤ M + 1 by the Cauchy-Schwarz inequality. From (3.14)  √    Z T  1 f M +1 M +1 1 √ E + P(τ ≤ T ) ≤ 1 + (3.15) + ¯ 1 ) 1τ 2 ≤t dt . f N q¯t (X N 0 t So far, we have not used the assumption q¯0 (x) = P(X0 = x) > 0 for any x ∈ LM . Since the jump intensity is bounded by

λf f ,

T − λf f

we necessarily have q¯t (x) ≥ e

q¯0 (x), for t ∈ [0, T ]. Thus,

there exists a constant C > 0 (depending on T and Linit ) such that q¯t 1(x) ≤ C for t ∈ [0, T ], x ∈ LM . Since τ 2 and τ 1 have the same law under P, this gives    √ Z T f M +1 M +1 1 1 √ P(τ ≤ t)dt , (3.16) +C P(τ ≤ T ) ≤ 1 + + f N N 0 and we easily conclude by Gronwall’s lemma.



Now, we can have an estimate of the accuracy given by the approximation (3.9). We assume that F : E → R is a bounded measurable function. Then, we know by the central limit theorem that ! N √ 1 X i i ¯t , Y¯t )t∈[0,T ] ) − E [F ((Xt , Yt )t∈[0,T ] )] → N (0, σ 2 ), N F ((X Q law N i=1

where σ 2 is the variance of F ((Xt , Yt )t∈[0,T ] ) under Q. Since F is bounded by a constant K > 0, we have by Proposition 10 N N N 1 X X K X 1 i,N i,N ¯ i , Y¯ i )t∈[0,T ] ) − Nα 1{τ i ≤T } , ≤ F ((X F ((X , Y ) ) t∈[0,T ] t t t t N N 1−α N i=1

i=1

i=1

which converges for the L1 -norm to 0 when N → +∞ for α < 1/2. Combining both results, we finally get a lower estimate of the convergence rate. 16

Corollary 11. Under the assumptions of Proposition 10, we have for any bounded measurable function F : E → R and any 0 < α < 1/2, N X i,N i,N α 1 N F ((Xt , Yt )t∈[0,T ] ) − EQ [F ((Xt , Yt )t∈[0,T ] )] → 0 in probability. N i=1

Remark 12. To prove Proposition 10, we have assumed that P(X0 = x) > 0 for any x ∈ LM . In fact, the same proof would work if we assumed that there existed x0 ∈ LM such that P(X0 ≥ x0 ) = 1 and P(X0 = x) > 0 for any x ≥ x0 . However, in practice, it would have been nice to treat the case P(X0 = x0 ) = 1 for some x0 ∈ LM , since we know at the beginning how many firms h have already i defaulted. Heuristically 1 ¯1 from (3.15), we may hope to have for large N that E q¯ (X¯ 1 ) 1τ 2 ≤t ≤ CP(τ 1 ≤ t) since X t t i h 1 2 1 and τ 2 are asymptotically independent, E q¯ (X ¯ 1 ) ≤ M + 1, and τ and τ have the same t t law. This would be enough to conclude. Unfortunately, despite our investigations, we have √ not been able to prove this formally. However, we still observe a convergence speed of 1/ N when X0 = 0 on our numerical experiments (Section 4).

4

Numerical results

In this section, we illustrate the theoretical results obtained in the previous sections. Let us recall that the local intensity (LI) model is a Markov chain with unit jumps occurring with the rate λ(t−, x) and that the SLI model is given by equation (1.3). First, we highlight that the LI and SLI models have the same marginal distributions but different laws as processes. Second, we study the convergence of the interacting particle system and obtain numerical simulations showing a central limit theorem. For our numerical experiments we consider two different models for the process Y described below and we assume that the local intensity λ and the function f are given by   ¯ 1− x λ(t, x) = λ M f (x) = (x ∨ f ) ∧ f where the number of defaultable entities M will be taken equal to M = 125 from now on. ¯ λ This choice of λ(t, x) corresponds to M independent default times with intensity M . We also assume that there is no default at the beginning, i.e. X0 = 0. 1. In the framework proposed by Lopatin and Misirpashaev (2008), the process Y is a continuous process satisfying a CIR type SDE p (4.1) dYt = κ(λ(t, Xt− ) − Yt ) + σ Yt dWt where κ, σ > 0. To sample such a process, we use the second-order discretization scheme for the CIR diffusion given in Alfonsi (2010).

2. In the framework of Arnsdorf and Halperin (2008), the process Y is no more continuous and may jump when a new default occurs. The process Y solves the following SDE dYt = −aYt log(Yt )dt + σYt dWt + γYt− dXt 17

(4.2)

where a, σ > 0 and γ ≥ 0. Remember that X only has positive jumps. For discretization purposes, note that between two jump times of X, log(Y ) solves the following Ornstein Uhlenbeck SDE 1 dZt = (−aZt + σ 2 )dt + σdWt . 2 Even though we could have sampled exactly in this particular case the Gaussian increments of Z, we discretize Y using the Euler scheme on Z in our simulation since we would make this choice for more general SDEs on Y . The LI model can be simulated very easily using a standard Monte–Carlo approach as it is sufficient to know how to sample a Poisson process with intensity λ; we do not need any interacting particle system.

4.1

Practical implementation

In this part, we describe our implementation of the particle system to sample from the distribution of (Xt , Yt )0≤t≤T . We recall that the process X has no continuous part and makes jumps of size 1. The process Y has a continuous part and may jump at the same times as X. T : sk = kh, 0 ≤ k ≤ D. Assume We consider a regular time grid of [0, T ] with step size h = D we have already discretized Y up to time sk , the discretization of Y at time sk+1 is built in the following way: • If the process X does not jump between time sk and time sk+1 we use the increment of a standard discretization scheme. • If the process X jumps at time s with sk < s < sk+1 , we proceed in three steps: apply the previous case between times sk and s, integrate the jump at time s and finally apply the previous case again between time s and sk+1 . This scheme ensures all the Y i are at least discretized on the regular grid {s0 , s1 , . . . , sD }. Computational complexity. Studying the complexity is of prime importance when proposing a numerical algorithm. On the one hand, there are D discretization times at which we recalculate the values of each Y i , which requires O(DN ) operations. On the other hand, the average total number of proposed jump dates (ie. the number of steps in loop 5) is given by the expectation of ¯¯ the underlying Poisson process at time T : N T λff . The computation cost of the ratio (4.3) ¯¯

is O(N ) and the complexity is O(N D + N 2 T λff ). Hence, for fixed model parameters, the overall complexity of our approach is bounded by O(DN + N 2 ). In practice, N is much larger than D, and the most computationally demanding part of the algorithm is the numerical approximation of the condition expectation involved in the jump intensity. The complexity of the interacting particle approach can be well improved if during the algorithm we keep track of the following two quantities involved in Equation (4.3) Dtj =

N X l=1

l )1{X l f (Yt−

t− =j}

and Ntj =

N X l=1

18

1{X l

t− =j}

j = 0, . . . , M.

(4.4)

Algorithm 1 (The naive particle system algorithm with complexity O(DN + N 2 ).) 1: t = 0 //Current date. 2: ti = 0 for all i = 1 . . . N //Last discretization date for the particle i. 3: Sample (X i , Y i ) independently according to the initial law. 4: sk = 0 // Last date on the regular grid: sk ≤ t < sk+1 . 5: while t ≤ T do ¯¯  6: t′ = t + E λff N //t′ is the potential next jump in the whole particle system. 7:

8: 9: 10: 11: 12: 13: 14: 15: 16:

17: 18: 19: 20: 21: 22: 23: 24: 25: 26:

while t′ > sk + h do sk = sk + h for i = 1 to N do Update the discretization of Y i from time ti to time sk . ti = s k end for t = sk end while I = uniform r.v. with values in {1, . . . , N }. //Index of the particle which may jump. Compute the conditional expectation PN l=1 1{X l =X I } . (4.3) E = PN l l=1 f (Y )1{X l =X I } f

R = λ¯ f¯λ(t′ , X I )f (Y I )E //Compute the acceptance ratio. U = uniform r.v. in [0, 1]. if U < R then //We accept the jump. Discretize Y I up to time t′ . Y I = Y I + γ(t′ , X I , Y I ) XI = XI + 1 tI = t′ end if t = t′ end while

Since the processes X l are unit-jump increasing these quantities clearly vanish for j > l . maxl Xt− Let t be the last proposed jump time of the particle system. These two vectors can be easily updated at time t′ which denotes the next possible jump time. If some ticks of the regular grid lie in [t, t′ ), we set t as the last discretization date in this interval and recompute vectors (Dtj )j and (Ntj )j using Equation (4.4). This happens D times in the algorithm and can be done with O(M + N ) operations as explained in Algorithm 2. Case 1: If the proposed jump at time t′ is not accepted, there is nothing to compute: Dtj′ = Dtj

and Ntj′ = Ntj .

19

Case 2: Otherwise, let I denote the index of the particle jumping  j j j+1 j+1 I I  Dt′ = Dt − f (Yt ); Dt′ = Dt + f (Yt′ ) Ntj′ = Ntj − 1; Ntj+1 = Ntj+1 + 1 ′   j Nt′ = Ntj ; Dtj′ = Dtj

at time t′ , we use if j = XtI , if j = XtI , otherwise.

(4.5)

Algorithm 2 (The improved particle system algorithm with complexity O(DN ).) 1: t = 0 2: ti = 0 for all i = 1 . . . N //Last discretization date for the particle i. 3: Sample (X i , Y i ) independently according to the initial law. 4: Set D l = N l = 0 for 0 ≤ l ≤ M . 5: for i = 1 to N do //Calculate D and N . We can directly set N 0 = N , D 0 = N f (Y0 ), l l D = N = 0 for 1 ≤ l ≤ M when the initial law is a Dirac mass at (0, Y0 ). i i i i 6: N X = N X + 1, D X = D X + f (Y i ) 7: end for 8: sk = 0 //Last date on the regular grid. 9: while t ≤ T do ¯¯  //t′ is the potential next jump in the whole particle system. 10: t′ = t + E λff N 11:

12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35:

while t′ > sk + h do sk = sk + h for i = 1 to N do Update the discretization of Y i from time ti to time sk . ti = s k end for t = sk Reinitialize D l = N l = 0 for 0 ≤ l ≤ M . for i = 1 to N do //Recalculate D and N. i i i i N X = N X + 1, D X = D X + f (Y i ) end for end while I = uniform r.v. with values in {1, . . . , N } //Index of the particle which may jump. f

XI

//Compute the acceptance ratio. R = λ¯ f¯λ(t′ , X I )f (Y I ) N X I D U = uniform r.v. in [0, 1] if U < R then //We accept the jump. I I I I N X = N X − 1, D X = D X − f (Y I ) Discretize Y I up to time t′ . Y I = Y I + γ(t′ , X I , Y I ) XI = XI + 1 I I I I N X = N X + 1, D X = D X + f (Y I ) tI = t′ end if t = t′ end while

20

One should notice that in cases 2 and 3, the updating cost does not depend on the size of the vectors (Dtj )j and (Ntj )j . Using these updating formulas, we can improve Algorithm 1 to obtain Algorithm 2. With this new algorithm we only have to compute D full approximations of the conditional expectation which has a unit a cost of O(M + N ) and the rest of the time ¯¯ we use the updating formulas (4.5), which happens on average less than N T λff times. Then, ¯¯

the overall cost of this new algorithm is O(D(M + N ) + N T λff ). For fixed model parameters, this complexity reduces to O(DN ). This new algorithm has a linear cost with respect to the number of particles. Thus, we managed to propose an interacting particle algorithm with the same cost as a crude Monte–Carlo method for SDEs since M is in practice fixed and much smaller than N . The CPU times of the two algorithms are compared in the following examples.

4.2

Marginal distributions

(a) Default distribution of the SLI model for Y given by Equation (4.2) with T = 1, Y0 = 1, a = 1, σ = ¯ = 2.5. 0.3, γ = 1, λ

(b) Default distribution of the SLI model for Y given by Equation (4.1) with T = 1, Y0 = 1, κ = 1, σ = ¯ = 2.5. 0.3, λ

(c) Default distribution of the LI model for T = ¯ = 2.5. 1, λ

(d) Binomial distribution function with parameters ¯ M , p = 1 − e−λT /M .

Figure 1: Comparison of the marginal distributions of the LI and SLI models. The simulations use N = 50, 000 samples, D = 100 and f = 1/3, f = 3.

21

In Figures 1, we draw the probability distribution of X for the LI and SLI models for both processes Y considered in this part. The comparison of these graphs highlights how the LI and SLI models effectively mimic their marginal distributions; their probability distributions look almost the same. Since the LI model corresponds to independent default times with ¯ λ intensity M , the default distribution at time T is actually the Binomial law with parameters   ¯

λ T M and p = 1 − exp − M

as we can see on Figure 1(d).

Model of Fig. 1(a) Model of Fig. 1(b)

Algorithm 1 239 229

Algorithm 2 1.51 1.49

Table 1: Comparison of the CPU times (in seconds) of the two algorithms with the same parameters as in Figure 1. We compare in Table 1 the computational times of the two algorithms and the gain obtained by the second approach is definitely outstanding. Algorithm 2 massively outperforms Algorithm 1 by a factor of 150. Of course, this gain will be all the more important as the number of particles increases. Given the impressive match of the marginal distributions, we would like to numerically investigate the difference between their distributions as processes. To do so, we have computed in each model the length of the longest interval during which X does not jump defined by τ = sup{t ∈ [0, T ] : ∃u ∈ [0, T − t], Xu+t− = Xu }

(4.6)

Note that with this definition, τ = T when X does not jump on the interval [0, T ]. The histograms of τ in the LI and SLI models are shown in Figure 2; we can see that, in the SLI, the length of the longest interval without jumps for X can be very small with a probability much higher than in the LI model (the l.h.s. of the histogram in the SLI model is fatter than in the LI model). This impression is reinforced by more quantitative observations. From the data used to plot these histograms, we have computed in Table 2 several values of the cumulative distribution function of the length of the longest interval without defaults both in the LI and SLI models. These quantities differ sufficiently to be numerically convinced that these two distributions do not match. SLI model LI model

P(τ ≤ T /4) 0.1911 (± 0.0033) 0.1645 (± 0.0033)

P(τ ≤ T /8) 0.0200 (± 0.0012) 0.0113 (± 0.0009)

Table 2: Some values of the distribution function of the length of the longest interval without ¯ = 2.5. The SLI model is defined as in Figure 2(a). The simulations jumps for T = 2 and λ use 50, 000 samples. The values between braces correspond to twice the standard deviation of the estimator.

4.3

Convergence of the interacting particle system

When introducing the interacting particle system, we emphasized that it was not only a theoretical tool but that it was also of practical interest as it satisfies a strong law of large numbers. From a numerical point of view, the efficiency of the particle system depends on the 22

(a) SLI model for Y given by Equation (4.2) with Y0 = 1, a = 1, σ = 0.3, γ = 3.

(b) LI model

Figure 2: Histogram of the length of the longest interval without X jumping for T = 2 and ¯ = 2.5. These histograms uses 50, 000 samples. λ

rate at which every particle converges to the invariant probability (see Section 3.5). In that section, we proved that this convergence rate was faster than N α for 0 < α < 1/2 where N is the number of particles. Now, we want to study this convergence rate in several examples: the first example only involves the marginal distribution of the particle system at maturity time, whereas the other two examples require the knowledge of the whole distribution and not only the marginal ones. Number of defaults distribution. First, we start with a simple example. We want to study the convergence of the estimator of P(XT = 3) computed on the particle system. We ran 5000 independent copies of the interacting particle systems and we computed the value of the estimator for each system. In Figure 3, we can see the centered and renormalised histogram of the values obtained for the empirical estimators. The histogram can be compared to the density of the standard normal distribution plotted as a solid line on Figure 3 and they match pretty well. This result suggests that a kind of central limit theorem should hold in practice even though we did not manage to prove such a result.

Asian option on the number of defaults. For our second example, we consider an Asian option on the default counting process whose price is given by   Z T 1 Xu du − K P =E T 0 + This price P will be approximated using the corresponding particle system estimate PN , where N is the number of particles. We are interested in the limiting distribution of PN . Because the process X has no continuous part and only makes jumps of size 1 and X0 = 0, the pathwise integral can be rewritten Z X 1 1 T t Xu du = XT − T 0 T t≤T s.t. Xt− 6=Xt

23

Figure 3: Centered and renormalised distribution of the estimation of P(XT = 3) using the particle system when Y is given by Equation (4.2) with T = 1, Y0 = 1, a = 1, σ = ¯ = 2.5, f = 1/3, f = 3 and N = 10000. The histogram was obtained using 5000 0.3, γ = 1, λ independent particle systems.

Hence, there is no need to approximate the integral, it can be computed exactly (up to the simulation of X). The example requires to sample the joint distribution of XT and the sum of the default times. On Figure 4, we have plotted the distribution of PN after renormalizing and centering. As before, the solid line is the standard normal density. We can see that the limiting distribution looks very much like the Gaussian distribution, but such an histogram does not enable to determine the rate of convergence to the limiting distribution. Actually, the rate of convergence p 1/2 . From a practical point of is given by the decrease rate of Var(PN ) = E(|PN − P |2 ) view we do not have access to P , so we have approximated it by the empirical mean Pˆ of the data set used to build the histogram of Figure 4. q We can see on Figure 5 that the rate of decrease of E(|PN − Pˆ |2 ) recalls the shape of a negative power function. Then, we have decided to compute the linear regression of − 21 log E(|PN − Pˆ |2 ) with respect to log(N ) on our simulations of PN for N varying from 100 to 10, 000 with a step size of 100, which gives a set of 100 data. The idea of the regression is to write 1 − log E(|PN − Pˆ |2 ) = α log(N ) + β + εN 2 P and to minimize the series N ε2N . The minimum is achieved for α = 0.5014, β = −0.1361 and the empirical variance of the sequence (εn ) is equal to √ 10−4 . This computation yields that the rate of convergence to the limiting distribution is N . It ensues from this result √ combined with the analysis of the histogram 4 that a Central Limit Theorem with rate N should hold. 24

Figure 4: Centered and renormalised distribution of the estimator PN of the Asian option price using the particle system when Y is given by Equation (4.1) with T = 2, Y0 = 1, a = 1, σ = 0.3, κ = 1, f = 1/3, f = 3 and N = 10000. The histogram was obtained using 10, 000 independent particle systems.

Figure 5: Convergence rate of

q E|PN − Pˆ |2 w.r.t the number of particles N .

25

Longest interval without jump. In this paragraph, we are interested in the convergence rate of the estimator of the length of the longest time interval with no jump. We recall the quantity of interest already defined in Equation (4.6) τ = sup{t ∈ [0, T ] : ∃u ∈ [0, T − t], Xu+t− = Xu } and we consider its particle system estimator τN . To numerically sample from the distribution of τ , we need to know the joint distribution of the jumping times of X, which is a prime case of a pathwise estimator. On Figure 6, we can see the centered and renormalized distribution of τN together with the standard normal density function plotted as a solid line. Again, the limiting distribution looks very much like a Gaussian distribution. Using these 5, 000 independent particle systems, each with 10, 000 particles, we can compute an approximation of τ , denote τˆ in the following. Now, we can run independent simulations of particle systems p with a number of particles N varying from 100 to 10, 000. We study the rate of decrease of E|τN − τˆ|2 , which p according to Figure 7 shows a negative power function shape. If we linearly regress log E|τN − τˆ|2 against log(N ) we find that we can write 1 − log E(|τN − τˆ|2 ) = α log(N ) + β + εN 2 The linear regression yields α = 0.4865, β = 1.016 and the empirical variance of the sequence (εn ) is equal to√0.0004. This regression yields that the rate of convergence to the limiting √ distribution is N , which focuses the existence of a Central Limit Theorem with rate N .

Figure 6: Centered and renormalised distribution of the estimator τˆ of τ using the particle system when Y is given by Equation (4.1) with T = 2, Y0 = 1, a = 1, σ = 0.3, κ = 1, f = 1/3, f = 3 and N = 10, 000. The histogram was obtained using 10, 000 independent particle systems.

26

Figure 7: Convergence rate of

5

p

E|τN − τˆ|2 w.r.t the number of particles N .

Conclusion

Local intensity models are wide spread for modelling the default counting process in credit risk. Recently, more sophisticated models with a stochastic factor involved in the intensity have been introduced in the literature. These stochastic local intensity models can be automatically calibrated to CDO tranche prices by properly choosing the local part of the intensity. This particular choice of the local intensity gives rise to a very specific family of SLI dynamics for which we have investigated the existence and uniqueness of solutions. This theoretical study has been carried out using particle systems, which turned out to be a clever tool for the numerical simulation of such dynamics. We have proved that particle Monte-Carlo algorithms based on this particle system approach almost surely converge. The theoretical study of the convergence rate enabled us to prove that the almost sure convergence took place at a rate faster that N α for any 0 < α < 1/2. Obtaining a Central Limit Theorem type result for such particle systems remains an open question, even though we could highlight such a behaviour in all our simulations. Last, we have shown that the interacting particle system can be sampled with a computational cost in O(N D), which is the same asymptotic cost as a Monte-Carlo algorithm for standard SDEs.

27

A

Proof of Proposition 3

The scheme of the proof is the following: • First, we prove that Ψ+ : t, x 7−→ Ψ(t, x+ ) is globally Lipschitz in x. Then, (E ′′ ) admits a unique solution on R+ . • Second, we prove that the solution satisfies ∀t ≥ 0, P (t) ≥ 0 and |P (t)| = 1. Step 1: Ψ+ is globally Lipschitz. Let us prove that for all (x, y) ∈ E × E there exists a constant K such that |Ψ+ (t, x) − Ψ+ (t, y)| ≤ K|x − y|. We have X X |(Ψ(t, x+ ))ij − (Ψ(t, y+ ))ij |. |Ψ+ (t, x) − Ψ+ (t, y)| = i∈LM j≥0

Bounding this quantity boils down to bound

X X

i∈LM j≥0

(A1 )ij := |

X k≥0

 (A1 )ij + 2f (j)λ(t, i)(A2 )ij where

µikj ((xik )+ − (yki )+ )|,

P∞ i P∞ i (xl )+ (yl )+ i i i l=0 l=0 (A2 )j := P∞ (xj )+ − P∞ (yj )+ . i i f (l)(x ) f (l)(y ) l=0

l +

l +

l=0

Bound for A1 . We easily get from Hypothesis 1 X X X XX X X |µikj ||(xik )+ − (yki )+ | = 2|µikk ||(xik )+ − (yki )+ | (A1 )ij ≤ i∈LM k≥0 j≥0

i∈LM j≥0

i∈LM k≥0

≤2

P∞

sup

i∈LM ,k∈N

|µikk ||x − y|.

(y i )

P∞

(y i )

l + l + (xij )+ and ± P∞ l=0 (xij )+ in Bound for A2 . To bound it, we introduce ± P∞ l=0 f (l)(xi ) f (l)(y i ) l +

l=0

l=0

l +

order to split (A2 )ij in three terms. Each of them is bounded in the following way

P∞ i P∞ i ∞ X (xij )+ P l=0 (xl )+ (xij )+ − P l=0 (yl )+ (xij )+ ≤ P |(xil )+ − (yli )+ |, ∞ ∞ ∞ f (l)(xi ) f (l)(xi ) f (l)(xi ) l=0

l +

l=0

l +

l + l=0 ∞ X (xi )+ P∞ j |(xil )+ i) f (l)(x l=0 l + l=0 l=0

P∞ i P∞ i P l=0 (yl )+ (xij )+ − P l=0 (yl )+ (xij )+ ≤ ∞ ∞ f (l)(xi ) i l=0 l=0 f (l)(yl )+ l + P∞ i P∞ i P l=0 (yl )+ (xij )+ − P l=0 (yl )+ (yji )+ ≤ ∞ ∞ f (l)(y i ) i l=0 l=0 f (l)(yl )+ l +

f f

− (yli )+ |,

1 i |(x )+ − (yji )+ |. f j

≤ λ(1 + 2 ff )|x − y|. Combining bounds on A1 and A2 ,   f i . we get |Ψ+ (t, x) − Ψ+ (t, y)| ≤ K|x − y|, where K = 2 sup |µkk | + 2λ 1 + 2 f i∈LM ,k∈N

Then,

P

i∈LM

P

i j≥0 f (j)λ(t, i)(A2 )j

Step 2: the solution is positive with norm 1. Let P (t) denote the unique solution of (E ′′ ). 28

First, we prove that ∀t ∈ R+ , ∀(i, j) ∈ LM × N, Pji (t) ≥ 0. Pji (t) satisfies (

(Pji )′ (t) =

P

λ(t,i−1)f (j) i−1 i i (t))+ − k6=j µkj (Pk (t))+ + 1{i≥1} ϕ((P (t))+ ,i−1) (Pj

Pji (0) = δix0 δjy0 .



λ(t,i)f (j)1{i≤M −1} ϕ((P (t))+ ,i)

 − µijj (Pji (t))+ ,

We assume that there exists t ≥ 0 such that Pji (t) < 0. We also introduce u := sup{s ≤ t : Pji (s) = 0}. Then, we integrate the above equation between u and t. We get (Pji )(t)

=

Z tX u k6=j

µikj (Pki (s))+ + 1{i≥1}

λ(s, i − 1)f (j) (P i−1 (s))+ ds. ϕ((P (s))+ , i − 1) j

The l.h.s. is strictly negative whereas the r.h.s. is non negative. Then, Pji (t) ≥ 0 for all t ∈ R+ . Second, ∀t ∈ R+ , P |P (t)| P = 1. Since P is non negative, (|P (t)|)′ = P Pwe prove i ′ i ′ i∈LM j≥0 (Pj ) (t). Moreover, i∈LM j≥0 (Pj ) (t) = 0. Then, |P (t)| = |P (0)| = 1.

B

Proof of Proposition 5

In fact, we can explicitly describe the solution of (3.3) as follows. We have (X0 , Y0 ) = (x0 , y0 ). Then, up to the first jump of X, Yt is necessarily defined as the unique strong solution of the following SDE Z t Z t σ(s, x0 , Ys )dWs , t ≤ T. b(s, x0 , Ys )ds + Yt = y 0 + 0

0

The process (Xt , t ≥ 0) only increases by jumps of size 1, and has at most M jumps since λ(t, M ) = 0. These jumps may only occur at the times T n and a jump do occur at time T n if the following condition is satisfied Un ≤

f λ(T n −, XT n − )f (YT n − ) . ϕQ (T n −, XT n − ) λf

Observe here that for any x ∈ LM , t 7→ ϕQ (t, x) is c`adl` ag since the process (X, Y ) is c`adl` ag under Q, and ϕQ (T n −, XT n − ) is well defined. If the jump occurs, we have YT n = YT n − + γ(T n −, XT n − , YT n − ), and Yt is defined up to the next jump of X as the unique strong solution of the SDE Z t Z t σ(s, XTn , Ys )dWs , Tn ≤ t ≤ T. b(s, XTn , Ys )ds + Yt = YT n + Tn

Tn

29

C

Proof of E[supt∈[0,T ] |Yt1,N |p ] < ∞

We assume that Y 1,N satisfies (3.6) and X 1,N is a Poisson process with intensity (3.5). Then Z t p Z t p  1,N p p p−1 1,N 1,N 1,N 1,N |Yt | ≤ 4 y0 + b(s, Xs , Ys )ds + σ(s, Xs , Ys )dWs 0 0 p  Z t 1,N 1,N + γ(s− , Xs1,N − , Ys− )dXs 0 Z u p  Z t p p p−1 1,N 1,N p−1 1,N 1,N b(s, Xs , Ys ) ds + sup ≤4 σ(s, Xs , Ys )dWs y0 + T u≤t

0

+ M p−1

Z

t

0

1,N p 1,N |γ(s− , Xs1,N − , Ys− )| dXs



0

,

since X 1,N jumps at most M times. The r.h.s being increasing w.r.t t, we can replace in the l.h.s. |Yt1,N |p by supu≤t |Yu1,N |p . Burkholder-Davis-Gundy inequality yields to p    Z t Z u 1,N 1,N p 1,N 1,N p/2−1 σ(s, Xs , Ys ) ds . σ(s, Xs , Ys )dWs ≤ Cp T E sup E u≤t

0

0

On the other hand, E

hR

t 1,N 1,N p 1,N − 0 |γ(s , Xs− , Ys− )| dXs

i



λf f

i h − , X 1,N , Y 1,N )|p ds E |γ(s − − 0 s s

Rt

since the jump intensity of X 1,N is upper bounded by λf f . Now since b, σ, and γ have a sub linear growth with respect to y (see Hypothesis (2)), we get  Z T 1,N p 1,N p 1 + E[sup |Yu | ]ds , E[sup |Yt | ] ≤ C t≤T

0

u≤s

which gives the result by Gronwall’s lemma.

D

Proof of Lemma 7

The proof of this lemma is done in two steps. First, we claim that (π N )N is tight if and only if the sequence ((Xt1,N , Yt1,N )t∈[0,T ] )N is tight. To check this, we first notice that P(E) is a Polish space. By Proposition 4.6 in M´el´eard (1996), (π N )N is tight if and only if (I(π N ))N is tight. Then Prohorov’s Theorem (tightness is equivalent to sequential compactness) gives with equation (3.7) the claim. Now, we must show that ((Xt1,N , Yt1,N )t∈[0,T ] )N is tight. We use Aldous’ criterion. First, we have to check that, for any ε > 0, there exists a constant K such that P(supt∈[0,T ] |Xt1,N | + |Yt1,N | > K) ≤ ε. This is trivial since Xt1,N is bounded and E[supt∈[0,T ] |Yt1,N |] < ∞ (see Appendix C). Second, we have to check that for any ε > 0, η > 0, there exist δ > 0 and n0 such that sup

sup

n≥n0 τ1 ,τ2 ∈T[0,T ] ,τ1 ≤τ2 ≤τ1 +δ

P(|(Xτ1,N , Yτ1,N ) − (Xτ1,N , Yτ1,N )| > η) < ε, 2 2 1 1

where T[0,T ] denotes the set of stopping times taking values in [0, T ]. We take |(x, y)| = max(|x|, |y|) and we assume without loss of generality that 0 < η < 1. For convenience, we 30

introduce νti,N =

P∞

k=1 1{T i,k ≤t} : this is a Poisson process 1,N ντ2 ≥ ντ1,N + 1 (X 1,N jumps in (τ1 , τ2 ] ) 1

the two cases: jump in (τ1 , τ2 ] ). We have:

λf f . We distinguish 1,N ντ1 (X 1,N does not

with intensity and ντ1,N = 2

 P(|(Xτ1,N , Yτ1,N ) − (Xτ1,N , Yτ1,N )| > η) ≤P ντ1,N ≥ ντ1,N +1 2 2 1 1 2 1 1,N Since P(ντ1,N > ντ1,N − ντ1,N 2 1 ) ≤ E[ντ2 1 ] =

 + P |Yτ1,N − Yτ1,N | > η, ντ1,N = ντ1,N := P1 + P2 . 2 1 2 1

λf f E[τ2

− τ1 ], we get P1 ≤ δ λf f .

Rτ Rτ P2 is bounded by P( τ12 b(s, Xs , Ys )ds + τ12 σ(s, Xs , Ys )dWs > η). Using Markov’s inequality, we get " Z 2 # Z τ2 τ2 1 1,N 1,N 1,N 1,N σ(s, Xs , Ys )dWs . b(s, Xs , Ys )ds + P2 ≤ 2 E η τ1 τ1 Moreover " Z E

2 # σ(s, Xs1,N , Ys1,N )dWs b(s, Xs1,N , Ys1,N )ds + τ1 τ1 " Z " Z 2 # 2 #! τ2 τ2 ≤ 2 E b(s, Xs1,N , Ys1,N )ds + E σ(s, Xs1,N , Ys1,N )dWs . τ2

Z

τ2

τ1

τ1

R 2 Rτ τ2 1,N 1,N On the one hand, we have τ1 b(s, Xs , Ys )ds ≤ δ τ12 C(1 + |Ys1,N |2 )ds for some constant inequality gives h R C > 0 by Hypothesis 2. iOn the other hR hand, Burkholder-Davis-Gundy’s i τ τ E | τ12 σ(s, Xs1,N , Ys1,N )dWs |2 ≤ CE τ12 1 + |Ys1,N |2 ds . Since E[supt≤T |Yt1,N |2 ] < ∞, we

get P2 ≤ ηC2 δ. Thus, combining the upper bounds on P1 and P2 , Aldous’ criterion is satisfied for δ := λf ε C . f

+

η2

References F. Abergel and R. Tachet. A nonlinear partial integro-differential equation from mathematical finance. Discrete Contin. Dyn. Syst., 27(3):907–917, 2010. doi: 10.3934/dcds.2010.27.907. C. Alexander and L. M. Nogueira. Stochastic local volatility. Proc. of the 2nd IASTED International Conference, pages 136–141, 2004. A. Alfonsi. High order discretization schemes for the CIR process: application to affine term structure and Heston models. Math. Comp., 79(269):209–237, 2010. doi: 10.1090/ S0025-5718-09-02252-2. A. Alfonsi. An introduction to the multiname modelling in credit risk. Recent Advancements in the Theory and Practice of Credit Derivatives, Eds T. Bielecki, D. Brigo, F. Patras, Bloomberg Press., 2011. M. Arnsdorf and I. Halperin. BSLP: Markovian bivariate spread-loss model for portfolio credit derivatives. J. Comput. Finance, 12(2):77–107, 2008. 31

R. Cont and A. Minca. Recovering portfolio default intensities implied by CDO quotes. Mathematical Finance, 23(1):94–121, 2013. doi: 10.1111/j.1467-9965.2011.00491.x. B. Dupire. Pricing with a smile. Risk, (7):18–20, 1994. J. Guyon and P. Henry-Labord`ere. The Smile Calibration Problem Solved. SSRN eLibrary, 2011. I. Gy¨ongy. Mimicking the one-dimensional marginal distributions of processes having an Itˆo differential. Probab. Theory Relat. Fields, 71(4):501–516, 1986. doi: 10.1007/BF00699039. J.-P. Lepeltier and B. Marchal. Probl`eme des martingales et ´equations diff´erentielles stochastiques associ´ees ` a un op´erateur int´egro-diff´erentiel. Ann. Inst. H. Poincar´e Sect. B (N.S.), 12(1):43–103, 1976. A. V. Lopatin and T. Misirpashaev. Two-dimensional Markovian model for dynamics of aggregate credit loss. In Econometrics and risk management, volume 22 of Adv. Econom., pages 243–274. Emerald/JAI, Bingley, 2008. doi: 10.1016/S0731-9053(08)22010-4. S. M´el´eard. Asymptotic behaviour of some interacting particle systems; McKean-Vlasov and Boltzmann models. In Probabilistic models for nonlinear partial differential equations (Montecatini Terme, 1995), volume 1627 of Lecture Notes in Math., pages 42–95. Springer, Berlin, 1996. doi: 10.1007/BFb0093177. V. Piterbarg. Markovian Projection Method for Volatility Calibration. SSRN eLibrary, 2006. P. Sch¨ onbucher. Portfolio losses and the term structure of loss transition rates: A new methodology for the pricing of portfolio credit derivatives. Working Paper, 2005. ´ ´ e de Probabilit´es de Saint-Flour A.-S. Sznitman. Topics in propagation of chaos. In Ecole d’Et´ XIX—1989, volume 1464 of Lecture Notes in Math., pages 165–251. Springer, Berlin, 1991. doi: 10.1007/BFb0085169.

32