2.1. Notation and preliminaries

The rate of going into default per obligor is $\lambda>0$ . This means that if there are still i obligors left (i.e. not in default), the time to the next default is exponentially distributed with mean $(\lambda i)^{-1}$ .

Recall that $p_n(u,t)$ is the probability of ruin before time t, starting with n obligors at time 0, given the initial reserve level is u. In our approach we (uniquely) characterize $p_n(u,t)$ through its double transform

\begin{equation*}\psi_n(\gamma) := \int_0^\infty \,{\mathrm{e}}^{-\gamma u} \int_0^\infty \vartheta \,{\mathrm{e}}^{-\vartheta t} p_n(u,t)\, {\mathrm{d}} t\, {\mathrm{d}} u=\int_0^\infty \,{\mathrm{e}}^{-\gamma u} p_n(u)\, {\mathrm{d}} u,\end{equation*}

where $p_n(u)$ can be interpreted as the probability of ruin before an exponentially distributed clock with mean $\vartheta^{-1}$ (which is sampled independently from anything else). The case of $t=\infty$ corresponds to $\vartheta\downarrow 0$ . The main result of this section is an expression (recursive in n) for $\psi_n(\gamma)$ : we express $\psi_n(\cdot)$ in terms of $\psi_{n-1}(\cdot)$ . Observe that we can equivalently write $p_n(u)$ as ${\mathbb P}(Z_n\geq u)$ , where $Z_n$ is the maximum of the net cumulative loss process (the net cumulative claim process, in the insurance context) over the above-mentioned exponentially distributed amount of time (i.e. with mean $\vartheta^{-1}$ ).

In practical settings one typically has $r>-\lambda\ell'(0)=\lambda\,{\mathbb E}L$ , so at any point in time ruin is rare, in the sense that the expected reserve increases as a function of time; to this end, note that when there are $i\in\{0,\ldots,n\}$ obligors left, the ‘local drift’ of the reserve process is $ri + \lambda i\,\ell'(0)>0$ .

2.2. Analysis

In this subsection we present a recursive scheme to evaluate $\psi_n(\gamma)$ . The main idea is to condition on the first event, being either the first default (which happens after an exponentially distributed time with mean $(\lambda n)^{-1})$ or the expiration of the exponential clock (which happens after an exponentially distributed time with mean $\vartheta^{-1})$ . If the former event happens to apply first, then we can still reach ruin, but now with $n-1$ obligors and an adapted initial reserve. If the latter events occur first, then we will not be facing ruin before the exponential clock expires. These ideas can be translated into mathematical terms as

(2) \begin{equation}p_n(u) = \int_0^\infty \lambda n\, {\mathrm{e}}^{- (\lambda n+\vartheta) t} \,{\mathbb P}(Z_{n-1}+L\geq u+rnt)\, {\mathrm{d}} t,\end{equation}

using the fact that the time to the first event is exponentially distributed with mean $(\lambda n+\vartheta)^{-1}$ , and that the first event is a default with probability $\lambda n/(\lambda n+\vartheta)$ .

We proceed by analyzing $\psi_n(\gamma)$ using relation (2), with the objective to express it in terms of $\psi_{n-1}(\cdot)$ . By a change of variable $v:=u+rnt$ , we obtain

\begin{align*}\psi_n(\gamma) &=\int_0^\infty \,{\mathrm{e}}^{-\gamma u} \int_0^\infty \lambda n\, {\mathrm{e}}^{- (\lambda n+\vartheta) t} \,{\mathbb P}(Z_{n-1}+L\geq u+rnt)\, {\mathrm{d}} t\, {\mathrm{d}} u\\[5pt] &=\dfrac{1}{rn} \int_0^\infty \,{\mathrm{e}}^{-\gamma u} \int_u^\infty \lambda n\, {\mathrm{e}}^{- (\lambda n+\vartheta) (v-u)/(rn)}\, {\mathbb P}(Z_{n-1}+L\geq v)\, {\mathrm{d}} v\, {\mathrm{d}} u.\end{align*}

The next step is to swap the order of the integrals, exploiting the fact that the integral over u allows an elementary solution:

\begin{align*}\dfrac{1}{rn}& \int_0^\infty \biggl(\int_0^v \,{\mathrm{e}}^{-\gamma u} \,{\mathrm{e}}^{ (\lambda n+\vartheta) u/(rn)} \, {\mathrm{d}} u\biggr) \lambda n\, {\mathrm{e}}^{- (\lambda n+\vartheta) v/(rn)}\, {\mathbb P}(Z_{n-1}+L\geq v)\, {\mathrm{d}} v\\[5pt] &=\dfrac{\lambda n}{\gamma rn -\lambda n-\vartheta}\int_0^\infty ({\mathrm{e}}^{- (\lambda n+\vartheta) v/(rn)}-{\mathrm{e}}^{-\gamma v})\, {\mathbb P}(Z_{n-1}+L\geq v)\, {\mathrm{d}} v.\end{align*}

In the last expression we see an object that resembles a Laplace transform, but observe that it features a complementary cumulative distribution function rather than a density. Recall, however, the standard identity

(3) \begin{equation}\int_0^\infty \,{\mathrm{e}}^{-\gamma u}\, {\mathbb P}(X\geq u) \,{\mathrm{d}} u =\dfrac{1}{\gamma}-\dfrac{1}{\gamma}\int_0^\infty \,{\mathrm{e}}^{-\gamma u}\, {\mathbb P}(X\in {\mathrm{d}} u)= \dfrac{1-{\mathbb E}\,{\mathrm{e}}^{-\gamma X}}{\gamma}.\end{equation}

In addition, using integration by parts, for the non-negative random variable $Z_{n-1}$ ,

(4) \begin{equation}{\mathbb E}\,{\mathrm{e}}^{-\gamma Z_{n-1}} = \int_0^\infty \,{\mathrm{e}}^{-\gamma x}\, {\mathbb P}(Z_{n-1}\in {\mathrm{d}} x) = 1-\gamma \int_0^\infty\, {\mathbb P}(Z_{n-1}>x)\,{\mathrm{e}}^{-\gamma x} \,{\mathrm{d}} x= 1-\gamma\psi_{n-1}(\gamma).\end{equation}

By the identity (3), and using the independence between the random variables $Z_{n-1}$ and L, we obtain, for any $\gamma\geq 0$ , with $d_n:=(\lambda n+\vartheta)/(rn)$ ,

\begin{align*}\psi_n(\gamma) &=\dfrac{\lambda n}{\gamma rn -\lambda n-\vartheta}\biggl(\dfrac{rn}{\lambda n+\vartheta}(1- {\mathbb E}\,{\mathrm{e}}^{-(\lambda n+\vartheta)/(rn)\,(Z_{n-1}+L)}) -\dfrac{1}{\gamma}(1- {\mathbb E}\,{\mathrm{e}}^{-\gamma\,(Z_{n-1}+L)})\biggr)\\[6pt] &=\dfrac{\lambda n}{\lambda n+\vartheta}\dfrac{1}{\gamma}+\dfrac{\lambda n}{\gamma rn -\lambda n-\vartheta}\biggl(\dfrac{{\mathbb E}\,{\mathrm{e}}^{-\gamma Z_{n-1}}\ell(\gamma)}{\gamma}-\dfrac{{\mathbb E}\,{\mathrm{e}}^{-d_n Z_{n-1}}\ell(d_n)}{d_n}\biggr),\end{align*}

which, by applying (4) and a few elementary algebraic steps, equals

\begin{equation*}\dfrac{\lambda n}{\lambda n+\vartheta}\dfrac{1}{\gamma}+\dfrac{\lambda n}{\lambda n +\vartheta-\gamma r n}\biggl(B\biggl(\dfrac{\lambda n +\vartheta}{rn},\psi_{n-1}\biggl(\dfrac{\lambda n +\vartheta}{rn}\biggr)\biggr)-B(\gamma,\psi_{n-1}(\gamma))\biggr),\end{equation*}

where we define

\begin{equation*}B(x,y):= \ell(x)\biggl(\dfrac{1}{x}-y\biggr).\end{equation*}

Note that we have expressed $\psi_n(\cdot)$ in terms of $\psi_{n-1}(\cdot)$ , so we would obtain a recursion if we had an explicit expression for $\psi_0(\cdot)$ . Recall that $\psi_0(\cdot)$ corresponds to ruin in the scenario without any obligor left. Obviously $p_0(u,t)\equiv 0$ for any u and t, entailing that $\psi_0(\gamma)\equiv 0$ for any value of $\gamma$ . It means that we can thus recursively compute $\psi_n(\gamma)$ . The theorem below summarizes the findings so far.

Theorem 1. For any $\gamma\geq 0$ and $n\in{\mathbb N}$ , we have the recursion

\begin{equation*}\psi_{n}(\gamma)= \dfrac{\lambda n}{\lambda n+\vartheta}\dfrac{1}{\gamma}+\dfrac{\lambda n}{\lambda n +\vartheta-\gamma r n}\biggl(B\biggl(\dfrac{\lambda n +\vartheta}{r n},\psi_{n-1}\biggl(\dfrac{\lambda n +\vartheta}{rn}\biggr)\biggr)-B(\gamma,\psi_{n-1}(\gamma))\biggr),\end{equation*}

where $\psi_0(\gamma)\equiv 0$ .

Remark 2. Upon inspecting the above proof, it is readily checked that we have not used the fact that the income rate is proportional to the number of obligors present; similarly, it is not crucial that the time to the next default when there are still i obligors is exponential with parameter $\lambda i$ . This effectively means that we can work with an income rate $r_i$ (rather than ri) and a default rate $\lambda_i$ (rather than $\lambda i$ ) during times that there are i obligors left. We thus obtain the recursion

\begin{equation*}\psi_{n}(\gamma)= \dfrac{\lambda_n}{\lambda_n+\vartheta}\dfrac{1}{\gamma}+\dfrac{\lambda_n}{\lambda_n +\vartheta-\gamma r_n}\biggl(B\biggl(\dfrac{\lambda_n +\vartheta}{r_n},\psi_{n-1}\biggl(\dfrac{\lambda_n +\vartheta}{r_n}\biggr)\biggr)-B(\gamma,\psi_{n-1}(\gamma))\biggr),\end{equation*}

where $\psi_0(\gamma)\equiv 0$ . We also remark that one can make the loss distribution dependent on the number of obligors in the system, by working with the transform $\beta_i(\cdot)$ when there are still i obligors that have not yet gone into default.

Remark 3. An interesting special case relates to the situation in which $r_n = r$ and $\lambda_n=\lambda$ , i.e. the conventional Cramér–Lundberg model. Sending $n\to\infty$ , one should recover the (transient version of the) Pollaczek–Khinchine formula. As an illustration, we show this computation for $\vartheta=0$ , writing a for $\lambda/r$ and assuming that $-a\ell'(0)<1$ . We obtain the following relation, with the limit of $\psi_n(\cdot)$ denoted by $\psi(\cdot)$ :

\begin{equation*}\psi(\gamma) = \dfrac{1}{\gamma}+\dfrac{a}{a-\gamma}(B(a,\psi(a))-B(\gamma,\psi(\gamma)).\end{equation*}

After some elementary algebra, it yields

\begin{equation*}1-\gamma\psi(\gamma) = \dfrac{\gamma}{\gamma-a+a\ell(\gamma)}\ell(a)(1-a\psi(a)).\end{equation*}

The constant $\ell(a)(1-a\psi(a))$ can be identified by observing that the left-hand side goes to 1 as $\gamma\downarrow 0$ ; hence an application of l’Hôpital’s rule yields that

\begin{equation*}\ell(a)(1-a\psi(a)) =\lim_{\gamma\downarrow 0}\dfrac{\gamma-a+a\ell(\gamma)}{\gamma} = 1+a\ell'(0).\end{equation*}

We conclude that

\begin{equation*}\psi(\gamma) =\dfrac{1}{\gamma}-\dfrac{1+a\ell'(0)}{\gamma-a+a\ell(\gamma)},\end{equation*}

which directly corresponds to the Pollaczek–Khinchine formula [Reference Asmussen and Albrecher5,Reference Dbicki and Mandjes11]. Our new results can thus be seen as a true generalization of the classical results from ruin theory.

Remark 4. The recursion featuring in Theorem 1 can be made more explicit when working with its generating function. To demonstrate this, we focus on the case of $\vartheta=0$ , $r_n=rn$ , and $\lambda_n=\lambda n$ . We have, again with $a=\lambda/r$ ,

\begin{equation*}\psi_n(\gamma) = \dfrac{1}{\gamma} +\dfrac{a}{a-\gamma}\biggl(\ell(a)\biggl(\dfrac{1}{a}-\psi_{n-1}(a)\biggr)-\ell(\gamma)\biggl(\dfrac{1}{\gamma}-\psi_{n-1}(\gamma)\biggr)\biggr).\end{equation*}

Thus, using $\psi_0(\gamma)=0$ , we obtain

\begin{align*}\Psi(z,\gamma)&:= \sum_{n=1}^\infty z^n\psi_n(\gamma)\\&=\sum_{n=1}^\infty z^n\dfrac{1}{\gamma}+z\dfrac{a}{a-\gamma} \sum_{n=1}^\infty z^{n-1}\biggl(\ell(a)\biggl(\dfrac{1}{a}-\psi_{n-1}(a)\biggr)-\ell(\gamma)\biggl(\dfrac{1}{\gamma}\psi_{n-1}(\gamma)\biggr)\biggr)\\&=\dfrac{z}{1-z}\dfrac{1}{\gamma}+z\dfrac{a}{a-\gamma}\biggl(\ell(a)\biggl(\dfrac{1}{a{(1-z)}}-\Psi(z,a)\biggr)-\ell(\gamma)\biggl(\dfrac{1}{\gamma{(1-z)}}-\Psi(z,\gamma)\biggr)\biggr) .\end{align*}

We conclude that

\begin{equation*}\Psi(z,\gamma)=\dfrac{1}{a-\gamma-za\,\ell(\gamma)}\,\biggl( \dfrac{z}{1-z} \dfrac{a-\gamma}{\gamma} + za\, \ell(a) \biggl(\dfrac{1}{a{(1-z)}}-\Psi(z,a)\biggr)-\dfrac{z}{1-z} \dfrac{ a\,\ell(\gamma)}{\gamma}\biggr).\end{equation*}

We are thus left with determining $\Psi(z,a)$ . For a and z fixed there is a unique positive $\gamma\equiv \gamma(z,a)$ for which the denominator equals 0 (as follows from the fact that $\nu(\gamma):=a-\gamma-za\,\ell(\gamma)$ is concave with $\nu(0)=a(1-z)>0$ and $\nu(\gamma)\to-\infty$ as $\gamma\to\infty$ ). We therefore have that in $\gamma\equiv \gamma(z,a)$ the numerator should equal 0 as well. This leads to

\begin{align*}\Psi(z,a) &= \dfrac{1}{a{(1-z)}}+\dfrac{1}{1-z}\dfrac{1}{\gamma(z,a)\,\ell(a)}\biggl(\dfrac{a-\gamma(z,a)}{a} -\ell(\gamma(z,a))\biggr)\\[6pt] &=\dfrac{1}{a{(1-z)}}+{\dfrac{a-\gamma(z,a)-a\,\ell(\gamma(z,a))}{(1-z)\,a\gamma(z,a)\,\ell(a)}}.\end{align*}

Combining the above, we have thus identified

\begin{equation*}\Psi(z,\gamma)=\dfrac{z}{1-z}\dfrac{1}{a-\gamma-za\,\ell(\gamma)}\biggl(\dfrac{a-\gamma -a\,\ell(\gamma)}{\gamma}-\dfrac{a-\gamma(z,a)-a\,\ell(\gamma(z,a))}{\gamma(z,a)}\biggr).\end{equation*}

Multiplying by $(1-z)$ , we obtain the transform at a geometrically distributed time with success probability z. Sending $z\uparrow 1$ , and realizing that $\gamma(1,a)=0$ , we recover the stationary result discussed in Remark 3.

3.1. Notation and preliminaries

Throughout this entire section we let the times-to-default $T_1,\ldots,T_n$ be non-negative i.i.d. random variables, with density $f(\cdot)$ and distribution function $F(\cdot)$ , distributed as a generic random variable T. Let $Z_n(t)$ be the net cumulative loss amount at time $t\geq 0$ , given that at time 0 there are $n\in{\mathbb N}$ obligors present. For $i=1,\ldots,n$ and $t\geq 0$ , we let $W_i(t)$ denote the net cumulative loss amount of the ith obligor at time t. By distinguishing between the scenario that obligor i has gone into default at time t and its complement, we can write $W_i(t)$ as

(5) \begin{equation}W_i(t):=1_{\{T_i\leq t\}}L_i-r\min\{T_i,t\}.\end{equation}

We define the moment generating function ${\mathbb E}\, {\mathrm{e}}^{\alpha L}$ of the loss L by $\bar\ell(\alpha):=\ell(-\alpha)$ . Then, due to fact that the obligors are statistically identical,

\begin{equation*}{\mathbb E}\,{\mathrm{e}}^{\alpha Z_n(t)} =({\mathbb E}\,{\mathrm{e}}^{\alpha W_1(t)} )^n.\end{equation*}

In addition, we can compute the moment generating function of the net loss amount of a single obligor at time t. By conditioning on the time-to-default, using (5),

\begin{align*}\omega_t(\alpha):={\mathbb E}\,{\mathrm{e}}^{\alpha W_1(t)}&=\bar\ell(\alpha)\int_0^t f(s) \,{\mathrm{e}}^{-r\alpha s} \,{\mathrm{d}} s +{\mathrm{e}}^{-r\alpha t}\int_t^\infty f(s) \,{\mathrm{d}} s\\[5pt] &=\bar\ell(\alpha)\int_0^t f(s) \,{\mathrm{e}}^{-r\alpha s} \,{\mathrm{d}} s +{\mathrm{e}}^{-r\alpha t} (1-F( t)).\end{align*}

For instance, in the special case when the times-to-default are exponentially distributed with mean $\lambda^{-1}$ , we have

\begin{equation*}\omega_t(\alpha)=(1-{\mathrm{e}}^{-(\lambda+r\alpha)t}) \dfrac{\lambda}{\lambda+r\alpha}\bar\ell(\alpha) + {\mathrm{e}}^{-(\lambda+r\alpha)t}.\end{equation*}

3.2. Large-deviations asymptotics

The goal of this subsection is to establish a limit theorem for our ruin probability, given that we start with n obligors and an initial capital reserve level $nu>0$ , as n grows large. In other words, we analyze how the probability

(6) \begin{equation}q_n(t):= p_n(nu,t)={\mathbb P}(\exists s\in[0,t]\colon Z_n(s) \geq nu)= {\mathbb P}\biggl(\exists s\in[0,t]\colon \sum_{i=1}^n W_i(s) \geq nu\biggr)\end{equation}

behaves as $n\to\infty$ . We do so under the evident ‘rarity condition’ that, for all $t\geq 0$ , ${\mathbb E}Z_n(t)$ is smaller than nu, or, equivalently,

\begin{equation*}\sup_{t\geq 0}({\mathbb P}(T\leq t) \,{\mathbb E}L - r\,{\mathbb E}\min \{T,t\} ) < u,\end{equation*}

where we use that ${\mathbb E}W_1(t)=\omega'_t(0) = {\mathbb P}(T\leq t) \,{\mathbb E}L - r\,{\mathbb E}\min \{T,t\}$ . We start by establishing a lower bound. The underlying principle is that the probability of a union of events is bounded from below by the probability of the most likely event among them. This entails that, for any $s\in[0,t]$ , we have $q_n(t)\geq \check q_n(s)$ , where

\begin{equation*}\check q_n(s):= {\mathbb P}\biggl(\sum_{i=1}^n W_i(s) \geq nu\biggr).\end{equation*}

Define the Legendre transform pertaining to $W_1(s)$ :

\begin{equation*}I(s):=\sup_{\alpha}(\alpha u - \log \omega_s(\alpha)).\end{equation*}

Because of the rarity condition $\omega'_s(0) <u$ for all $s\geq 0$ , we can restrict ourselves to maximizing over $\alpha>0$ only; we define $\alpha^\star(s):= \arg\sup_{\alpha}(\alpha u - \log \omega_s(\alpha))$ . By Cramér’s theorem [Reference Dembo and Zeitouni13], we immediately have that, for any $s\in[0,t]$ ,

(7) \begin{equation}\liminf_{n\to\infty}\dfrac{1}{n}\log q_n(t) \geq\liminf_{n\to\infty}\dfrac{1}{n}\log \check q_n(s) =-\alpha^\star(s) u + \log \omega_s(\alpha^\star(s))=-I(s).\end{equation}

We also define

\begin{equation*}t^\star:= \arg\inf_{s\in[0,t]}I(s),\end{equation*}

which has the informal interpretation of the most likely time $Z_n(\cdot)$ exceeds nu. From the fact that the lower bound (7) applies for any $s\in[0,t]$ , we thus obtain that

\begin{equation*}\liminf_{n\to\infty}\dfrac{1}{n}\log q_n(t) \geq-\inf_{s\in[0,t]}I(s) = -I(t^\star).\end{equation*}

We proceed by proving that $-I(t^\star)$ is also an upper bound on the decay rate of $q_n(t)$ . The first step is to realize that ruin occurs at the default time of one of the n obligors. As a consequence, we can rewrite $q_n$ in terms of the union of n events:

\begin{equation*}q_n(t)={\mathbb P}\biggl(\exists j\in\{1,\ldots,n\}\colon T_j\in[0,t], \sum_{i=1}^n W_i(T_j) \geq nu\biggr),\end{equation*}

instead of the union of uncountably many events featuring in the representation (6). By the union bound, we obtain that this probability $q_n(t)$ is majorized by $n\hat q_n(t)$ , where

\begin{equation*}\hat q_n(t):={\mathbb P}\biggl(T_1\in[0,t], \sum_{i=1}^n W_i(T_1) \geq nu\biggr).\end{equation*}

As $n^{-1} \log n \to 0$ , it suffices to prove that $\limsup_{n\to\infty} n^{-1} \log \hat q_n(t)\leq -I(t^\star)$ . To this end, by conditioning on $T_1$ ,

\begin{equation*}\hat q_n(t) = \int_0^t f(s) \,{\mathbb P}\biggl( \sum_{i=2}^{n} W_i(s) +L_1 -rs \geq nu\biggr) \,{\mathrm{d}} s.\end{equation*}

Then observe that the $W_i(T_1)$ are dependent, but once conditioned on $T_1=s$ they have become independent. The next step is to apply the Markov inequality: for any $\alpha\geq 0$ , with $L_1$ being independent from $W_2(s),\ldots,W_{n}(s)$ ,

\begin{align*}{\mathbb P}\biggl( \sum_{i=2}^{n} W_i(s) +L_1 -rs \geq nu \biggr)& = {\mathbb P}\biggl( \exp\biggl({\alpha \sum_{i=2}^{n} W_i(s) +\alpha L_1}\biggr) \geq \exp({\alpha (nu+rs)})\biggr) \\&\leq (w_s(\alpha))^{n-1}\bar\ell(\alpha) \,{\mathrm{e}}^{-\alpha(nu+rs)}\\&\leq (w_s(\alpha))^{n-1}\bar\ell(\alpha) \,{\mathrm{e}}^{-\alpha(n-1)u}.\end{align*}

Upon combining the above, we have thus found that for any $\alpha(\cdot)\geq 0$ ,

\begin{equation*}\limsup_{n\to\infty}\dfrac{1}{n}\log \hat q_n(t) \leq \limsup_{n\to\infty}\dfrac{1}{n}\log\int_0^t f(s) \,(w_s(\alpha(s)))^{n-1} \bar\ell(\alpha(s))\,{\mathrm{e}}^{-\alpha(s)\,(n-1)u} \,{\mathrm{d}} s.\end{equation*}

Recall that, for any $t\geq 0$ , $I(t)=\alpha^\star(t) u - \log \omega_t(\alpha^\star(t))$ . Plugging in $\alpha(\cdot)=\alpha^\star(\cdot),$ we thus obtain, in the second inequality using the definition of $t^\star$ ,

(8) \begin{align}\notag \limsup_{n\to\infty}\dfrac{1}{n}\log \hat q_n (t)&\leq \limsup_{n\to\infty}\dfrac{1}{n}\log\int_0^t f(s)\, \bar\ell(\alpha^\star(s))\,{\mathrm{e}}^{-(n-1)I(s)} \,{\mathrm{d}} s\nonumber\\ &\leq \limsup_{n\to\infty}\dfrac{1}{n}\log\int_0^t f(s)\, \bar\ell(\alpha^\star(s))\,{\mathrm{e}}^{-(n-1)I(t^\star)} \,{\mathrm{d}} s\nonumber\\ &=-I(t^\star)+ \limsup_{n\to\infty}\dfrac{1}{n}\log\int_0^t f(s)\, \bar\ell(\alpha^\star(s))\, {\mathrm{d}} s.\end{align}

Observe that we are done if we succeed in proving that the second term in (8) is 0, for which it suffices to prove that the integral appearing in this term is finite. To this end, first observe that, with $\tau(\alpha):= {\mathbb E}\,{\mathrm{e}}^{\alpha T}$ ,

\begin{equation*}\lim_{t\to\infty} \omega_t(\alpha) = \bar\ell(\alpha) \tau(-r\alpha)=:\Delta(\alpha),\end{equation*}

so that $\alpha^\star(\infty)$ solves $\Delta'(\alpha)/\Delta(\alpha) = u$ .

Under Assumption 1 we have $\sup_{s\in[ 0,t]}\alpha^\star(s)\leq M$ for some finite M. Note that this holds whenever the function $\alpha^\star(\cdot)$ is continuous, whereas for $t=\infty$ we additionally require $\alpha^\star(\infty)<\infty$ . With this assumption in place and using that $\alpha\mapsto \bar\ell(\alpha)$ is increasing, we conclude that

\begin{equation*}\int_0^t f(s)\, \bar\ell(\alpha^\star(s))\, {\mathrm{d}} s \leq \bar\ell(M) \int_0^t f(s)\, {\mathrm{d}} s \leq \bar\ell(M)<\infty.\end{equation*}

Summarizing, we have shown

\begin{equation*}\limsup_{n\to\infty}\dfrac{1}{n}\log q_n(t) \leq -I(t^\star).\end{equation*}

We have arrived at the following result.

Theorem 2. As $n\to\infty$ , under Assumption 1,

\begin{equation*}\dfrac{1}{n}\log q_n(t) \to -I(t^\star).\end{equation*}

3.3. Efficient simulation

As the above theorem only provides us with the logarithmic asymptotics of $q_n$ , it is inherently imprecise. For instance, if the true asymptotic shape of $q_n$ is $n^\alpha \,\exp({-nI(t^\star)})$ for some $\alpha\in{\mathbb R}$ , or $\exp(n^\eta) \exp({-nI(t^\star)})$ for some $\eta\in (0,1)$ , the effect of the $\alpha$ and $\eta$ is not visible. One can get accurate estimates in an efficient way, however, applying importance sampling. Below we present an importance sampling algorithm, which we prove to be logarithmically efficient.

The key idea is that we decompose our rare-event probability $q_n$ into n rare-event probabilities, which we will be dealing with separately. We write

(9) \begin{equation}q_n(t) = \sum_{j=1}^n q_{nj}(t),\end{equation}

where

\begin{equation*}q_{nj}(t):={\mathbb P}(\mathscr{F}_j),\:\:\:\:\mathscr{F}_j:=\mathscr{E}_j\cap \bigcap_{i=1}^{j-1} \mathscr{E}_i^{\rm c},\:\:\:\:\mathscr{E}_j:= \biggl\{T_j\in[0,t], \sum_{i\not=j} W_i(T_j) +L_j -rT_j \geq nu\biggr\};\end{equation*}

the validity of (9) is due to the events $\mathscr{F}_j$ being (by construction) disjoint, while the union of the $\mathscr{E}_j$ equals the union of the $\mathscr{F}_j$ . The problem of efficiently estimating $q_n(t)$ thus reduces to the problem of efficiently estimating each of the $q_{nj}(t)$ (and adding up all the resulting estimates).

Fix a $j\in\{1,\ldots,n\}$ and focus on the estimation of $q_{nj}$ . We now define an importance sampling probability measure ${\mathbb Q}$ .

• Under ${\mathbb Q}$ the density of $T_j$ remains $f(\cdot)$ .

• Conditionally on $T_j=s$ , the moment generating function of $L_j$ becomes

  \begin{equation*}\bar\ell^{\mathbb Q}(\alpha) = \dfrac{\bar\ell(\alpha+\alpha^\star(s)}{\bar\ell(\alpha^\star(s))}.\end{equation*}
  Sampling $L_j$ from ${\mathbb Q}$ amounts to sampling from an exponentially twisted version of the actual distribution. This is a standard procedure in applied probability; for many frequently used distributions the twisted distribution remains within the same class of distributions, but with different parameters. For instance, the $\alpha$ -twisted version of an exponentially distributed random variable with parameter $\mu$ corresponds to an exponentially distributed random variable with parameter $\mu-\alpha$ (requiring that $\alpha\in[0,\mu)$ ).

• Conditionally on $T_j=s$ , the moment generating function of $W_i(s)$ (for $i\not= j$ ) becomes

  (10) \begin{equation}\omega_s^{\mathbb Q}(\alpha):= \dfrac{\omega_s(\alpha+\alpha^\star(s))}{\omega_s(\alpha^\star(s))}.\end{equation}
  To decide whether the event $\mathscr{F}_j$ applies, we have to sample the default times $T_i$ and (if $T_i<t$ ) the losses $L_i$ , for $i\not=j$ , in accordance with (10). This can be done as follows. By (10), the exponentially twisted version of $W_i(s)$ has the moment generating function
  \begin{align*} \omega_s^{\mathbb Q}(\alpha) = \dfrac{1}{\omega_s(\alpha^\star(s))}\biggl(\int_0^s f(v)\, {\mathrm{e}}^{-(\alpha+\alpha^\star(s))\, rv} \bar\ell(\alpha+\alpha^\star(s)) \,{\mathrm{d}} v + \int_s^\infty f(v)\, {\mathrm{e}}^{-(\alpha+\alpha^\star(s))\, rs} \,{\mathrm{d}} v\biggr).\end{align*}
  From this identity we observe that the $L_i$ can be sampled from a distribution with moment generating function $\bar\ell^{\mathbb Q}(\cdot)$ , as defined above, whereas the density $f^{\mathbb Q}(\cdot)$ of the $T_i$ (for $i\not=j$ ) becomes
  \begin{equation*}f^{\mathbb Q}(v) =\dfrac{1}{\omega_s(\alpha^\star(s))}f(v)( {\mathrm{e}}^{-\alpha^\star(s)\, rv} \bar\ell(\alpha^\star(s))1_{\{v\leq s\}} + {\mathrm{e}}^{-\alpha^\star(s)\, rs} 1_{\{v>s\}}).\end{equation*}

We proceed by detailing the importance-sampling-based simulation procedure, and establishing its asymptotic efficiency. To this end, we first observe that a generic sample of the likelihood ratio, say $\mathscr{L}_j$ , has the form

\begin{equation*}{{\mathrm{e}}^{-\alpha^\star(T_j)\,L_j}}\cdot{\bar\ell(\alpha^\star(T_j))}\prod_{i\not=j} ({{\mathrm{e}}^{-\alpha^\star(T_j)\,W_i(T_j)}}\cdot{\omega_{T_j}(\alpha^\star(T_j))}).\end{equation*}

Recall that on the event $\mathscr{F}_j$ we have $ \sum_{i\not= j} W_i(T_j) +L_j -rT_j \geq nu$ . As a consequence, on the event $\mathscr{F}_j$ the likelihood ratio $\mathscr{L}_j$ is majorized by

\begin{align*} &{\mathrm{e}}^{-\alpha^\star(T_j)(nu + r T_j)}\cdot{\bar\ell(\alpha^\star(T_j))}\cdot(\omega_{T_j}(\alpha^\star(T_j)))^{n-1}\\&\quad \leq {\mathrm{e}}^{-\alpha^\star(T_j)(n-1)u}\cdot{\bar\ell(\alpha^\star(T_j))}\cdot(\omega_{T_j}(\alpha^\star(T_j)))^{n-1}\\ &\quad =\bar\ell(\alpha^\star(T_j))\,{\mathrm{e}}^{-(n-1)\,I({T_j})} \\ &\quad \leq \bar\ell(M)\,{\mathrm{e}}^{-(n-1)\,I({T_j})}\\&\quad \leq \bar\ell(M)\,{\mathrm{e}}^{-(n-1)\,I({t^\star})},\end{align*}

with M as defined in Section 3.2 (where we let Assumption 1 be in force). We thus find that, with $\mathscr{I}_j$ denoting the indicator function of $\mathscr{F}_j$ , the almost sure inequality $\mathscr{L}_j\,\mathscr{I}_j \leq \bar\ell(M)\,{\mathrm{e}}^{-(n-1)\,I({t^\star})}$ , and therefore

\begin{equation*}\sum_{j=1}^n \mathscr{L}_j\,\mathscr{I}_j \leq n \,\bar\ell(M)\,{\mathrm{e}}^{-(n-1)\,I({t^\star})}.\end{equation*}

Evidently, to obtain an estimator with good precision, we have to repeat the above experiment sufficiently often. Suppose, for each $j\in\{1,\ldots,n\}$ , we perform $N\in{\mathbb N}$ independent trials. The corresponding likelihood ratios are denoted by $\mathscr{L}_{j,k}$ , and the indicator functions are $\mathscr{I}_{j,k}$ , with $j\in\{1,\ldots,n\}$ and $k\in\{1,\ldots,N\}$ . Our estimator thus becomes

\begin{equation*}\xi_N:=\dfrac{1}{N} \sum_{k=1}^N \sum_{j=1}^n \mathscr{L}_{j,k}\,\mathscr{I}_{j,k},\end{equation*}

which is (by construction) unbiased. The next step is to analyze the performance of this estimator. To this end, we observe in relation to its second moment that

\begin{equation*}{\mathbb E}_{\mathbb Q}\biggl(\biggl(\sum_{j=1}^n \mathscr{L}_{j}\,\mathscr{I}_{j}\biggr)^2\biggr) \leq n^2 \,(\bar\ell(M))^2\,{\mathrm{e}}^{-2(n-1)\,I({t^\star})},\end{equation*}

with ${\mathbb E}_{\mathbb Q}(\cdot)$ denoting expectation under ${\mathbb Q}$ . We find the following upper bound for the second moment:

\begin{equation*}\limsup_{n\to\infty}\dfrac{1}{n}\log {\mathbb E}_{\mathbb Q}\biggl(\biggl(\sum_{j=1}^n \mathscr{L}_{j}\,\mathscr{I}_{j}\biggr)^2\biggr) \leq -2 I(t^\star).\end{equation*}

By Theorem 2, and using the fact that variances are non-negative, we also have the corresponding lower bound:

\begin{align*} \liminf_{n\to\infty}\dfrac{1}{n}\log {\mathbb E}_{\mathbb Q}\biggl(\biggl(\sum_{j=1}^n \mathscr{L}_{j}\,\mathscr{I}_{j}\biggr)^2\biggr)& \geq \liminf_{n\to\infty}\dfrac{2}{n}\log {\mathbb E}_{\mathbb Q}\biggl(\sum_{j=1}^n \mathscr{L}_{j}\,\mathscr{I}_{j}\biggr) \\&= \liminf_{n\to\infty}\dfrac{2}{n}\log q_n(t) \\ &= -2 I(t^\star).\end{align*}

The above bounds lead to the following conclusion, which in practical terms entails that the number of runs needed to obtain an estimate of a given relative precision grows sub-exponentially in n. For the definition of logarithmic efficiency, and related performance notions in rare-event simulation, we refer to [6, Ch. VI].

3.4. Uniform bound

Intrinsic drawbacks of the large-deviations asymptotics are that they only kick in for large n, and they provide us with the decay rate only. This motivates the search for a uniform upper bound on the ruin probability $p_n(u,t)$ . The result is a Lundberg-type inequality derived along the same lines as was done in [5, Section XIII.5a] for the conventional Cramér–Lundberg model in which claims (or losses in the credit context) arrive according to a fixed-intensity Poisson process. We focus on the situation that when there are n obligors the time to the first default is exponentially distributed with mean $\lambda_n^{-1}$ and the income rate is $r_n$ . Let $\gamma_n$ be the positive solution for $\gamma$ in

\begin{equation*}\bar\ell(\gamma) \dfrac{\lambda_n}{\lambda_n+\gamma r_n}=1.\end{equation*}

Theorem 4. Suppose that $\gamma_n$ is non-increasing in n. Then

\begin{equation*}p_n(u,t)\leq p_n(u,\infty)\leq {\mathrm{e}}^{-\gamma_n u}.\end{equation*}

Proof. It is evident that $p_n(u,t)\leq p_n(u,\infty)$ . Let $Y_n$ be distributed as $L-r_n\,T_1$ , where $T_1$ is assumed exponentially distributed with mean $\lambda_n^{-1}$ (independent of L). Conditioning on $Y_n$ immediately yields

\begin{equation*}p_n(u,\infty)=\mathbb{P}(Y_{n}>u)+\int_{-\infty}^u p_{n-1}(u-y,\infty)\,{\mathbb P}({Y_{n}}\in {\mathrm{d}} y).\end{equation*}

We claim that this implies $p_n(u,\infty)\leq {\mathrm{e}}^{-\gamma_n u}$ . The proof is by induction. First note that the claim holds true for $n=0$ as $p_0(u,\infty)=0$ for all $u> 0$ . Assuming the inequality holds true for $n-1$ ,

\begin{align*}p_{n}(u,\infty)&\leq\mathbb{P}(Y_{n}>u)+\int_{-\infty}^u \,{\mathrm{e}}^{-\gamma_{n-1} (u-y)}\,{\mathbb P}({Y_{n}}\in {\mathrm{d}} y)\\&\leq\mathbb{P}(Y_{n}>u)+\int_{-\infty}^u \,{\mathrm{e}}^{-\gamma_{n} (u-y)}\,{\mathbb P}({Y_{n}}\in {\mathrm{d}} y)\\&\leq {\mathrm{e}}^{-\gamma_{n} u}\int^{\infty}_u \,{\mathrm{e}}^{\gamma_{n} y}\,{\mathbb P}({Y_{n}}\in {\mathrm{d}} y)+\int_{-\infty}^u \,{\mathrm{e}}^{-\gamma_{n} (u-x)}\,{\mathbb P}({Y_{n}}\in {\mathrm{d}} y)\\&={\mathrm{e}}^{-\gamma_{n} u}\,{\mathbb E}\,{\mathrm{e}}^{\gamma_{n} Y_n}\\&={\mathrm{e}}^{-\gamma_{n} u}\,\bar\ell(\gamma_n) \dfrac{\lambda_n}{\lambda_n+\gamma_n r_n}\\&={\mathrm{e}}^{-\gamma_n u},\end{align*}

where in the second inequality we have used the fact that $\gamma_n$ is non-increasing in n. □

4.1. Non-default losses

In this subsection we consider the following extension of the model analyzed in Section 2 (or indeed the more general one featured in Remark 2). Next to losses due to defaults (happening at a Poisson rate $\lambda_n$ with the losses having Laplace transform $\ell(\cdot)$ when n obligors are present), there are losses that do not correspond to defaults (happening at a Poisson rate $\lambda^\circ_n$ with the losses having Laplace transform $\ell^\circ(\cdot)$ when n obligors are present).

We again start our derivations by conditioning on the first event, being the first default, the first loss (not leading to default), or the expiration of the exponential clock. If a default happens first, then we can still reach ruin, but now with $n-1$ obligors and an adapted initial reserve. If the first event is a loss that does not correspond to a default, then we can still reach ruin with n obligors but an adapted initial reserve. If the exponential clock expires, then we will not be facing ruin.

This idea can be formalized as follows. With $L^\circ$ denoting a generic random variable corresponding to a non-default loss, we obtain the relation

\begin{align*}p_n(u) &= \int_0^\infty \,{\mathrm{e}}^{- (\bar \lambda_n+\vartheta) t}(\lambda_n\, {\mathbb P}(Z_{n-1}+L\geq u+r_nt)+ \lambda^\circ_n\, {\mathbb P}(Z_n+L^\circ\geq u+r_nt)) \,{\mathrm{d}} t.\end{align*}

Going through the same type of computations as those relied on in Section 2, we end up with a relation between $\psi_n(\cdot)$ and $\psi_{n-1}(\cdot)$ . More specifically, for any $\gamma\geq 0$ , using the notation $\bar{\lambda}_n=\lambda_n+\lambda_n^\circ$ , we find that

(11) \begin{align}\notag \psi_{n}(\gamma)&= \dfrac{\bar\lambda_n}{\bar\lambda_n+\vartheta}\dfrac{1}{\gamma} + \dfrac{\lambda_n}{\bar\lambda_n+\vartheta-\gamma r_n}\biggl(B\biggl(\dfrac{\bar\lambda_n+\vartheta}{r_n},\psi_{n-1}\biggl(\dfrac{\bar\lambda_n+\vartheta}{r_n}\biggr)\biggr)-B(\gamma,\psi_{n-1}(\gamma))\biggr)\\[6pt] &\quad\, + \dfrac{\lambda_n^\circ}{\bar\lambda_n+\vartheta-\gamma r_n}\biggl(B^\circ\biggl(\dfrac{\bar\lambda_n+\vartheta}{r_n},\psi_{n}\biggl(\dfrac{\bar\lambda_n+\vartheta}{r_n}\biggr)\biggr)-B^\circ(\gamma,\psi_{n}(\gamma))\biggr),\end{align}

where $B^\circ(\cdot\,,\cdot)$ is defined as $B(\cdot\,,\cdot)$ but with $\ell(\cdot)$ replaced by $\ell^\circ(\cdot)$ . Unfortunately this relation between $\psi_n(\cdot)$ and $\psi_{n-1}(\cdot)$ cannot be directly written in terms of an explicit recursion (as opposed to the model without non-default losses; see Theorem 1). The $\psi_n(\cdot)$ , however, can still be found recursively, using the following procedure.

To this end, we start by defining the (yet unknown) constants

\begin{equation*}A_n:=B^\circ\biggl(\dfrac{\bar\lambda_n+\vartheta}{r_n},\psi_{n}\biggl(\dfrac{\bar\lambda_n+\vartheta}{r_n}\biggr)\biggr).\end{equation*}

Then, using that $\psi_0(\cdot)\equiv 0$ , observe that $\psi_1(\gamma)$ obeys

(12) \begin{align}\notag \psi_1(\gamma)& = \dfrac{\bar\lambda_1}{\bar\lambda_1+\vartheta}\dfrac{1}{\gamma}+\dfrac{\lambda_1}{\bar\lambda_1+\vartheta-\gamma r_1}\biggl(\dfrac{r_1}{\bar\lambda_1+\vartheta}\ell\biggl(\dfrac{\bar\lambda_1+\vartheta}{r_1}\biggr)-\dfrac{\ell(\gamma)}{\gamma}\biggr)\nonumber\\[5pt] &\quad\, + \dfrac{\lambda_1^\circ}{\bar\lambda_1+\vartheta-\gamma r_1}\biggl(B^\circ\biggl(\dfrac{\bar\lambda_1+\vartheta}{r_1},\psi_{1}\biggl(\dfrac{\bar\lambda_1+\vartheta}{r_1}\biggr)\biggr)-B^\circ(\gamma,\psi_{1}(\gamma))\biggr).\end{align}

We can rewrite (12), for a known function $F(\cdot)$ , as

\begin{equation*}\psi_1(\gamma) = F(\gamma) +\dfrac{\lambda_1^\circ}{\bar\lambda_1+\vartheta-\gamma r_1}\biggl(A_1-\ell^\circ(\gamma)\biggl(\dfrac{1}{\gamma}-\psi_1(\gamma)\biggr)\biggr),\end{equation*}

which can be rearranged to

\begin{equation*}1-\gamma \psi_1(\gamma) = 1 -\dfrac{\gamma F(\gamma)(\bar\lambda_1+\vartheta-\gamma r_1)+\gamma\lambda_1^\circ A_1-\lambda_1^\circ\ell^\circ(\gamma)}{\bar\lambda_1-\lambda_1^\circ\ell^\circ(\gamma)+\vartheta-\gamma r_1}.\end{equation*}

As we know that $1-\gamma \psi_1(\gamma)$ is a Laplace transform, its value should be between 0 and 1 for any $\gamma\geq 0$ . Hence any zero of the denominator is necessarily also a zero of the numerator. It is standard to verify that the numerator has a single positive zero, say $\bar\gamma$ . Then it follows that

\begin{equation*}A_1 = \dfrac{\ell^\circ(\bar\gamma)}{\bar{\gamma}}-F(\bar\gamma)\dfrac{\bar\lambda_1+\vartheta-\bar\gamma r_1}{\lambda_1^\circ}.\end{equation*}

Now that we have found $A_1$ and hence $\psi_1(\gamma)$ , we can identify $A_2$ and $\psi_2(\gamma)$ along the same lines: we first express $\psi_2(\gamma)$ in terms of $A_2$ using (11), and then identify $A_2$ using the fact that the zero of the denominator (which we know to equal $\bar\lambda_2 -\lambda_2^\circ\ell^\circ(\gamma)+\vartheta-\gamma r_2$ ) is a zero of the numerator as well. Continuing this procedure, all $\psi_n(\gamma)$ (and constants $A_n$ ) can be found.

4.2. Markov modulation

In the models discussed so far the individual obligors are independent. In reality they may be affected by common external factors, to be thought of as the ‘state of the economy’, and hence behave dependently. In this subsection we consider a model in which a particular dependence structure is incorporated, via the mechanism of Markov modulation (also known as regime switching).

We start by describing the model. Let $(J(t))_{t\geq 0}$ be an irreducible continuous-time Markov process living on $\{1,\ldots,d\}$ . We let $q_{jk}\geq 0$ (for $j\not=k$ ) denote the transition rate from state j to state k, and $q_j:=-q_{jj}=\sum_{k\not=j} q_{jk}$ . Let $r_{nj}$ be the rate at which the surplus process increases when there are n obligors and the background process is in state j, let $\lambda_{nj}$ be the corresponding hazard rate of the time to the next default, and let $\ell_j(\cdot)$ be the Laplace transform of the loss (with the associated generic random variable being denoted by $L_j$ ).

Let $T_n$ be the minimum of the time of the first default and the expiration of an exponential clock of rate $\vartheta$ . Let

\begin{equation*}R(T_n):=\int_0^{T_n} r_{nJ(t)} \,{\mathrm{d}} t\end{equation*}

denote the increase of the surplus process until $T_n$ . We start by analyzing the distribution of $R(T_n)$ through the object

\begin{equation*}F_{i,j,n}(x):= {\mathbb P}_i(R(T_n)\geq x, J(T_n)=j) := {\mathbb P}(R(T_n)\geq x, J(T_n)=j\,|\, J(0)=i).\end{equation*}

Using the standard ‘Markovian reasoning’, i.e. by distinguishing between all possible events in a (small) time interval of length $\Delta$ and using the memory-less property, we obtain the relation, as $\Delta\downarrow 0$ ,

\begin{equation*}F_{i,j,n}(x) = \sum_{k\not = j}F_{i,k,n}(x)\,q_{kj}\Delta + F_{i,j,n}(x-r_j\Delta)(1- (q_j+\lambda_{nj}+\vartheta))+o(\Delta).\end{equation*}

Subsequently subtracting $F_{i,j,n}(x-r_j\Delta)$ from both sides, dividing by $\Delta$ , and taking the limit $\Delta\downarrow 0$ , we end up with a system of linear differential equations:

\begin{equation*}F'_{i,j,n}(x) = \sum_{k=1}^d F_{i,k,n}(x)\,q_{kj} + F_{i,j,n}(x)\,(\lambda_{nj}+\vartheta).\end{equation*}

For given i and n, this is a system of d coupled linear differential equations, that can be solved in the standard manner; the resulting structure depends on the multiplicities of the eigenvalues. Henceforth we assume that its solution is such that the corresponding density obeys

\begin{equation*}{\mathbb P}_i(R(T_n)\in {\mathrm{d}} x, J(T_n)=j) = \sum_{k=1}^d \xi_{i,j,k,n} \,{\mathrm{e}}^{-\zeta_{k,n}x},\end{equation*}

but a similar analysis can be done if the terms in the right-hand side of the previous display also involve polynomial factors (as a consequence of the multiplicities of some of the eigenvalues being larger than one).

The key observation is the identity

\begin{align*}{\mathbb P}_i(Z_n\geq u) &= \sum_{j=1}^d \dfrac{\lambda_{nj}}{\lambda_{nj}+\vartheta}\int_0^\infty\, {\mathbb P}_j(Z_{n-1}\in {\mathrm{d}} z) \, {\mathbb P}_i(L_j\geq R(T_n)+u-z, J(T_n)=j)\\&= \sum_{j=1}^d \dfrac{\lambda_{nj}}{\lambda_{nj}+\vartheta}\int_0^\infty\int_0^\infty\, {\mathbb P}_j(Z_{n-1}\in {\mathrm{d}} z)\, {\mathbb P}(L_j\geq x+u-z) \sum_{k=1}^d \xi_{i,j,k,n} \,{\mathrm{e}}^{-\zeta_{k,n}x}\, {\mathrm{d}} x\\&= \sum_{j=1}^d \dfrac{\lambda_{nj}}{\lambda_{nj}+\vartheta}\int_0^\infty\, {\mathbb P}_j(Z_{n-1}+L_j\geq x+u)\sum_{k=1}^d \xi_{i,j,k,n} \,{\mathrm{e}}^{-\zeta_{k,n}x}\, {\mathrm{d}} x {.}\end{align*}

Therefore, using the now familiar steps concerning a change of variables and swapping the order of integration,

\begin{align*}\psi_{ni}(\gamma)&:=\int_0^\infty \,{\mathrm{e}}^{-\gamma u}\, {\mathbb P}_i(Z_n\geq u)\, {\mathrm{d}} u\\&=\sum_{j=1}^d \dfrac{\lambda_{nj}}{\lambda_{nj}+\vartheta}\int_0^\infty\int_0^\infty \,{\mathrm{e}}^{-\gamma u}\, {\mathbb P}_j(Z_{n-1}+L_j\geq x+u)\sum_{k=1}^d \xi_{i,j,k,n} \,{\mathrm{e}}^{-\zeta_{k,n}x}\, {\mathrm{d}} x\, {\mathrm{d}} u\\&=\sum_{j=1}^d \dfrac{\lambda_{nj}}{\lambda_{nj}+\vartheta}\int_0^\infty\int_u^\infty \,{\mathrm{e}}^{-\gamma u}\, {\mathbb P}_j(Z_{n-1}+L_j\geq v)\sum_{k=1}^d \xi_{i,j,k,n} \,{\mathrm{e}}^{-\zeta_{k,n}(v-u)}\, {\mathrm{d}} v\, {\mathrm{d}} u\\&=\sum_{j=1}^d \dfrac{\lambda_{nj}}{\lambda_{nj}+\vartheta}\int_0^\infty\sum_{k=1}^d \xi_{i,j,k,n}\biggl(\int_0^v \,{\mathrm{e}}^{-\gamma u} \,{\mathrm{e}}^{\zeta_{k,n}u}\, {\mathrm{d}} u\biggr)\,{\mathbb P}_j(Z_{n-1}+L_j\geq v) \,{\mathrm{e}}^{-\zeta_{k,n}v}\, {\mathrm{d}} v\\&=\sum_{j=1}^d \dfrac{\lambda_{nj}}{\lambda_{nj}+\vartheta}\int_0^\infty\sum_{k=1}^d \xi_{i,j,k,n}\dfrac{{\mathrm{e}}^{-\zeta_{k,n}v}-{\mathrm{e}}^{-\gamma v}}{\gamma- \zeta_{k,n}}\,{\mathbb P}_j(Z_{n-1}+L_j\geq v)\, {\mathrm{d}} v.\end{align*}

From now on we can follow the approach presented in Section 2: the last expression in the previous display can be expressed in terms of $\psi_{n-1,j}(\cdot)$ , for $j=1,\ldots,d$ . We thus end up with a vector-valued recursion. As the derivation is fully analogous to that corresponding to the non-modulated case, we omit the details.

4.3. Brownian perturbations

We proceed by making the model more realistic, by allowing the process to evolve, between defaults, as Brownian motion rather than a deterministic drift. The parameters of this Brownian motion depend on the number of obligors that have not yet gone into default, say with drift coefficient $r_i$ and variance coefficient $\sigma_i^2$ when there are i obligors left. In this section the time between the ith and $(i+1)$ th default is exponentially distributed with mean $\lambda_i^{-1}$ .

Considering a Brownian motion with parameters r and $\sigma^2$ over an interval with exponentially distributed length with mean $\lambda^{-1}$ , we know the following from Wiener–Hopf theory.

• The maximum value $M^+$ achieved is exponentially distributed with the parameter

  \begin{equation*}\nu^+\equiv \nu^+(r,\sigma^2,\lambda):= \dfrac{\sqrt{r^2+2\lambda\sigma^2}}{\sigma^2}-\dfrac{r}{\sigma^2}.\end{equation*}

• The (absolute value of the) amount by which the process goes down after the maximum is achieved until the end of the exponentially distributed interval, say $M^-$ , is exponentially distributed with the parameter

  \begin{equation*}\nu^-\equiv \nu^-(r,\sigma^2,\lambda):= \dfrac{\sqrt{r^2+2\lambda\sigma^2}}{\sigma^2}+\dfrac{r}{\sigma^2}.\end{equation*}

• The random variables $M^+$ and $M^-$ are independent. The rates $\nu^+$ and $\nu^-$ are the roots of the equation $\lambda+r\alpha-\tfrac{1}{2}\alpha^2\sigma^2=0$ .

Now define $\nu_n^\pm:= \nu^\pm(-r_n,\sigma_n^2,\lambda_n+\vartheta)$ ; note that the first parameter is $-r_n$ rather than $r_n$ , as we consider the event of the cumulative claim process exceeding the value u (i.e. the reserve level dropping below 0). As before, we set up a relation between $\psi_n(\cdot)$ and $\psi_{n-1}(\cdot)$ . Realize that, due to the Brownian term, ruin can occur before the exponential clock (with parameter $\vartheta$ ) expires; this happens with probability $ {\mathrm{e}}^{-\nu_n^+ u}$ . Following the approach we have been using in the case without the Brownian term, we thus obtain the relation

\begin{equation*}p_n(u) = {\mathrm{e}}^{-\nu_n^+ u} + I_n(u,\vartheta),\end{equation*}

where

\begin{align*}I_n(u,\vartheta)&:=\int_0^u \int_0^\infty\nu_n^+ \,{\mathrm{e}}^{-\nu_n^+ v}\nu_n^- \,{\mathrm{e}}^{-\nu_n^- w}\dfrac{\lambda_n}{\lambda_n+\vartheta}\,{\mathbb P}(Z_{n-1}+L\geq u-v+w)\, {\mathrm{d}} w\, {\mathrm{d}} v\\[5pt] &=\dfrac{\lambda_n}{\lambda_n+\vartheta}\int_0^u \int_{u-v}^\infty\nu_n^+ \,{\mathrm{e}}^{-\nu_n^+ v}\nu_n^- \,{\mathrm{e}}^{-\nu_n^- (z-u+v)}\,{\mathbb P}(Z_{n-1}+L\geq z)\, {\mathrm{d}} z\, {\mathrm{d}} v.\end{align*}

The next step is to evaluate $\psi_n(\gamma)$ , by multiplying $p_n(u)$ by ${\mathrm{e}}^{-\gamma u}$ and integrating over $u\in[0,\infty)$ . Interchanging the order of the integrals so that the ‘easy’ integration (i.e. over u) can be done first, we obtain

\begin{align*}\int_0^\infty &{\mathrm{e}}^{-\gamma u} I_n(u,\vartheta) \,{\mathrm{d}} u \\[5pt] &= \dfrac{\lambda_n}{\lambda_n+\vartheta}\int_0^\infty\int_0^\infty\int_v^{z+v} \,{\mathrm{e}}^{-\gamma u} \,\nu_n^+ \,{\mathrm{e}}^{-\nu_n^+ v}\nu_n^- \,{\mathrm{e}}^{-\nu_n^- (z-u+v)}\,{\mathbb P}(Z_{n-1}+L\geq z)\, {\mathrm{d}} u\, {\mathrm{d}} v\, {\mathrm{d}} z\\[5pt] &= \dfrac{\lambda_n}{\lambda_n+\vartheta}\int_0^\infty\int_0^\infty {\nu_n^-}\,{\mathrm{e}}^{-\gamma v}\dfrac{{\mathrm{e}}^{-\nu_n^- z}-{\mathrm{e}}^{-\gamma z}}{\gamma-\nu_n^-}\nu_n^+\,{\mathrm{e}}^{-\nu_n^+ v}\,{\mathbb P}(Z_{n-1}+L\geq z)\, {\mathrm{d}} v\, {\mathrm{d}} z\\[5pt] &= \dfrac{\lambda_n}{\lambda_n+\vartheta} \dfrac{\nu_n^-\nu_n^+}{(\gamma-\nu_n^-)(\gamma+\nu_n^+)} \int_0^\infty({\mathrm{e}}^{-\nu_n^- z}-{\mathrm{e}}^{-\gamma z}) \,{\mathbb P}(Z_{n-1}+L\geq z)\, {\mathrm{d}} z.\end{align*}

Performing the same steps as in the proof of Theorem 2.2, as before relying on the identities (3) and (4) in combination with the independence of L and $Z_{n-1}$ , after some standard algebra we obtain the following result.

Theorem 5. For any $\gamma\geq 0$ and $n\in{\mathbb N}$ , we have the recursion

\begin{align*}\psi_n(\gamma)&=\dfrac{1}{\gamma+\nu_n^+}+\dfrac{\lambda_n}{\lambda_n+\vartheta}\dfrac{1}{\gamma+\nu_n^+}{\dfrac{\nu_n^+}{\gamma}}\\[5pt] &\quad\, -\dfrac{\lambda_n}{\lambda_n+\vartheta} \dfrac{\nu_n^-\nu_n^+}{(\gamma-\nu_n^-)(\gamma+\nu_n^+)}( B(\nu_n^-,\psi_{n-1}(\nu_n^-))-B(\gamma,\psi_{n-1}(\gamma))),\end{align*}

where $\psi_0(\gamma)\equiv 0$ .

Remark 6. In Theorem 5 we can simplify

\begin{equation*}\dfrac{\lambda_n}{\lambda_n+\vartheta} \dfrac{\nu_n^-\nu_n^+}{(\gamma-\nu_n^-)(\gamma+\nu_n^+)}=\dfrac{\lambda_n}{\lambda_n+\vartheta+r_n\gamma-\frac{1}{2}\gamma^2\sigma_n^2},\end{equation*}

using that $\nu_n^+$ and $\nu_n^-$ solve $(\lambda_n+\vartheta)+r_n\alpha-\tfrac{1}{2}\alpha^2\sigma_n^2=0$ .

4.4. Multiple groups

To make the model more realistic, one could work with multiple (heterogeneous) groups of obligors. Suppose there are $G\in\mathbb{N}$ groups of obligors with initially $n_j$ obligors in group $j\in\{1,\ldots,G\}$ ; write ${\boldsymbol n}=(n_1,\ldots, n_G)$ . We consider the multi-group counterpart of the base model of Section 2: each obligor in group j has a time-to-default that is exponentially distributed with rate $\lambda_j$ . The losses at default per obligor in group j are i.i.d. random variables with Laplace transform $\ell_j(\cdot)$ ; in addition these per-group sequences are assumed independent. The income per unit time for this group is $r_j i$ when there are $i\in\{1,\ldots, n_j\}$ obligors that have not yet gone into default.

The company’s capital reserve is given by the sum of the reserves of the individual groups; its initial level is $u>0$ . Let $\psi_{{\boldsymbol n}}(\gamma)$ denote the double transform of the probability of ruin over an exponentially distributed interval (with, as usual, mean $\vartheta^{-1}$ ), given that there are $n^j$ obligors in group j that have not yet gone into default. Then, by the same argumentation as before, we find, for ${\boldsymbol n}$ component-wise at least equal to 1, and with ${\boldsymbol e}_j$ the jth unit vector

\begin{align*}\psi_{\boldsymbol n}(\gamma)&=\sum_{j=1}^G\dfrac{\lambda_j n_j}{\sum_{k=1}^G\lambda_k n_k+\vartheta}\dfrac{1}{\gamma}+\sum_{j=1}^G\dfrac{\lambda_j n_j}{\sum_{k=1}^G\lambda_k n_k +\vartheta-\gamma r_j n_j}\\[6pt] &\quad\, \times\biggl(B_j\biggl(\dfrac{\lambda_j +\vartheta/n_j}{r_j},\psi_{{\boldsymbol n}-{\boldsymbol e}_j}\biggl(\dfrac{\lambda_j +\vartheta/n_j}{r_j}\biggr)\biggr)-B_j{(\gamma,\psi_{{\boldsymbol n}-{\boldsymbol e}_j}(\gamma))}\biggr),\end{align*}

where

\begin{equation*}B_j(x,y):= \ell_j(x)\biggl(\dfrac{1}{x}-y\biggr).\end{equation*}

We have thus expressed $\psi_{\boldsymbol n}(\gamma)$ as a linear function of $\psi_{{\boldsymbol n}-{\boldsymbol e}_1}(\gamma)$ up to $\psi_{{\boldsymbol n}-{\boldsymbol e}_G}(\gamma)$ . A similar recursive relation be found if some of the entries of ${\boldsymbol n}$ equal 0. Given that $\psi_{{\boldsymbol 0}}(\gamma)=0$ , with ${\boldsymbol 0}$ denoting the G-dimensional all-zeros vector, we have thus devised a procedure to identify $\psi_{\boldsymbol n}(\gamma)$ .

5.1. Exponentially distributed losses

In Section 2.2 the focus was on finding an expression for the double transform $\psi_n(\gamma)$ , which can then be inverted numerically. In Section 3 we presented a couple of other approaches: asymptotics, an efficient importance sampling algorithm, and bounds. In this section we present an alternative technique, namely an iterative procedure that directly provides the ruin probabilities $p_n(u,t)$ themselves. We consider the model variant in which the default rate and the income rate are $\lambda_i$ and $r_i$ , respectively, during time periods in which there are i obligors left.

As in Section 2.2, the idea is to condition on the first default. We thus obtain, with $W(\cdot)$ as introduced in Section 3, the following recursive relation:

(13) \begin{align}\notag p_n(u,t)&=\int_0^t \lambda_n \,{\mathrm{e}}^{- \lambda_n s}\, {\mathbb P}\biggl(\sup_{0\leq v\leq t-s} \sum_{i=1}^{n-1}W_i(v)+L\geq u+r_ns\biggr)\, {\mathrm{d}} s\\ \notag&=\int_0^t \lambda_n \,{\mathrm{e}}^{- \lambda_ns} \, {\mathrm{d}} s-\int_0^t \lambda_n \,{\mathrm{e}}^{- \lambda_n s}\, {\mathbb P}\biggl(\sup_{0\leq v\leq t-s} \sum_{i=1}^{n-1}W_i(v)+L\leq u+r_ns\biggr)\, {\mathrm{d}} s\\&=1-{\mathrm{e}}^{-\lambda_n t}-\int_0^t\int_0^{u+r_n s} \lambda_n \,{\mathrm{e}}^{- \lambda_n s}(1-p_{n-1}(u+r_n s-x,t-s))\,\mathbb{P}(L\in {\mathrm{d}} x)\, {\mathrm{d}} s. \end{align}

When there is only one obligor left, there is only one scenario leading to ruin: default should take place before the exponential clock (with mean $\vartheta^{-1}$ ) expires and the loss should be sufficiently large. In other words,

\begin{align}\notag p_1(u,t)&= \int_0^t\int^\infty_{u+r_1s} \lambda_1 \,{\mathrm{e}}^{- \lambda_1 s}\, \mathbb{P}(L\in {\mathrm{d}} x)\, {\mathrm{d}} s= \int_0^t \lambda_1 \,{\mathrm{e}}^{- \lambda_1 s}\, \mathbb{P}(L\geq u+r_1s)\, {\mathrm{d}} s\end{align}

From this point on we focus on the case of exponentially distributed claims with mean $\mu^{-1}$ , i.e. ${\mathbb P}(L\geq x)={\mathrm{e}}^{-\mu x}$ . We readily obtain

\begin{equation*}p_1(u,t)=\int_0^t \lambda_1 \,{\mathrm{e}}^{- \lambda_1 s} \,{\mathrm{e}}^{-\mu (u+r_1s)}\, {\mathrm{d}} s=\dfrac{\lambda_1 \,{\mathrm{e}}^{-\mu u}}{\lambda_1+\mu r_1}(1-{\mathrm{e}}^{-(\lambda_1+\mu r_1) t}).\end{equation*}

We can thus obtain $p_2(u,t)$ applying numerical integration to (13) with $n=2$ . Continuing along these lines, $p_n(u,t)$ can be numerically evaluated for higher values of n.

We now point out how to evaluate the large-deviations asymptotics that were presented in Section 3.2, in the case of exponentially distributed claims. The moment generating function of $W_1(s)$ for $\alpha<\mu$ is given by

\begin{equation*}\omega_s(\alpha)=(1-{\mathrm{e}}^{-(\lambda+ r\alpha)s}) \dfrac{\lambda}{\lambda+r\alpha}\dfrac{\mu}{\mu-\alpha} + {\mathrm{e}}^{-(\lambda+r\alpha)s},\end{equation*}

whereas for $\alpha\geq \mu$ the moment generating function is infinite. We continue by computing the mean net loss corresponding to a single obligor (as a function of time):

\begin{align*}m(s)&:={\mathbb E}W_1(s) \\ &= \dfrac{1}{\mu}(1-{\mathrm{e}}^{-\lambda s}) - r\int_0^s u\, \lambda \,{\mathrm{e}}^{-\lambda v} \,{\mathrm{d}} v- rs\int_s^\infty \lambda \,{\mathrm{e}}^{-\lambda v} \,{\mathrm{d}} v\\&=\biggl(\dfrac{1}{\mu}-\dfrac{r}{\lambda}\biggr)(1-{\mathrm{e}}^{-\lambda s}).\end{align*}

Henceforth we will assume $u>m(\infty)$ , or equivalently $\lambda-r\mu<\lambda\mu u$ , to make sure the event under consideration is rare.

The Legendre transform pertaining to $W_1(s)$ reads

\begin{equation*}I(s):=\sup_{0<\alpha<\mu}(\alpha u - \log \omega_s(\alpha));\end{equation*}

we can rule out $\alpha\geq \mu$ as $\omega_s(\alpha)=\infty$ for these $\alpha$ . Because the first-order condition does not allow an explicit solution, one cannot write I(s) in closed form. Two boundary cases can be dealt with explicitly, though. It is first observed that, letting $\omega'_{s,1}(\alpha)$ denote the derivative of $\omega_s(\alpha)$ with respect to $\alpha$ , and letting $\omega'_{s,2}(\alpha)$ be the derivative of $\omega_s(\alpha)$ with respect to s,

(14) \begin{align} I'(s)& = \dfrac {\mathrm{d}} { {\mathrm{d}} s} (\alpha^\star(s) u - \log \omega_s(\alpha^\star(s)))\notag \\[5pt] &=\dfrac{ {\mathrm{d}} \alpha^\star(s)}{ {\mathrm{d}} s}\biggl(u- \dfrac{\omega'_{s,1}(\alpha^\star(s))}{\omega_s(\alpha^\star(s))}\biggr)- \dfrac{\omega'_{s,2}(\alpha^\star(s))}{\omega_s(\alpha^\star(s))}\notag \\[5pt] &=- \dfrac{\omega'_{s,2}(\alpha^\star(s))}{\omega_s(\alpha^\star(s))},\end{align}

where the last equality is due to the definition of $\alpha^\star(s)$ . By an elementary computation,

(15) \begin{equation}\omega'_{s,2}(\alpha)= \biggl(\dfrac{\lambda\mu}{\mu-\alpha}-(\lambda+r\alpha)\biggr)\,{\mathrm{e}}^{-(\lambda +r\alpha)s}=\dfrac{r\alpha^2+\lambda\alpha -r\mu\alpha}{\mu-\alpha}\,{\mathrm{e}}^{-(\lambda +r\alpha)s}.\end{equation}

We observe that the Legendre transform I(s) is decreasing in s whenever $\alpha^*(s)>\mu-{\lambda}/{r}$ .

• For $s=0$ , we immediately see that $\omega_0(\alpha) = 1$ for all $\alpha$ , so $\alpha^\star(0)=\mu$ and $I(0) = \mu u$ . In addition, we obtain by some straightforward algebra that

  \begin{equation*}I'(0) = - \lim_{\alpha\uparrow \mu} \dfrac{\omega'_{0,2}(\alpha)}{\omega_0(\alpha)}=-\infty.\end{equation*}

• For $s=\infty$ ,

  \begin{equation*}I(s) = \sup_{0<\alpha<\mu} \kappa(\alpha),\:\:\:\:\kappa(\alpha):=\alpha u - \log (\lambda\mu) +\log(\lambda+r\alpha)+\log(\mu-\alpha).\end{equation*}
  Observe that $\kappa(\cdot)$ is concave, with $\kappa'(0)>0$ (under the assumption $u>m(\infty)$ ) and $\kappa(\alpha)\to-\infty$ as $\alpha\uparrow\mu$ . In other words, $\kappa(\cdot)$ attains a maximum in $(0,\mu)$ . The first order condition, determining $\alpha^\star(\infty)$ , is
  \begin{equation*}u=\dfrac{1}{\mu-\alpha}-\dfrac{r}{\lambda+r\alpha},\end{equation*}
  or equivalently
  \begin{equation*}ru\alpha^2 +((\lambda-r\mu )u+2r)\alpha - \lambda\mu(u-m(\infty))=0.\end{equation*}
  As $\lambda\mu (u-m(\infty))>0$ , this equation has a positive and negative root. Consequently $\alpha^\star(\infty)$ is the positive root, that is,
  \begin{equation*}\alpha^\star(\infty) = \dfrac{-2r-\lambda u+r\mu u +\sqrt{4r^2+\lambda^2 u^2+2r\lambda\mu u^2 + r^2\mu^2 u^2}}{2ru},\end{equation*}
  so $I(\infty) = \kappa(\alpha^\star(\infty))$ . Next, we want to find the sign of I(s) as $s\to\infty$ . Based on (14) and (15), this is the sign of $-r\alpha^\star(\infty)-\lambda+r\mu$ . Using the explicit solution of $\alpha^\star(\infty)$ , it requires some straightforward calculus to verify that this leads to a negative sign, i.e. I(s) is decreasing as $s\to\infty$ , if and only if $ \lambda- r\mu>-\lambda\mu u$ .

5.2. Numerical example

For the numerical results we have used a set-up that aligns with that considered in [Reference Asmussen3].

• We consider the case when both the income rates $r_i$ and the default intensity $\lambda_i$ are linear in the number of obligors i that have not yet gone into default. We let the proportionality constants be $r= 1$ and $\lambda=0.9$ , respectively. In other words, when there are i obligors in the system that have not yet gone into default, the income rate is given by i and the default intensity rate by $0.9\,i$ .

• The losses are exponentially distributed with parameter $\mu=1$ .

With these parameter settings the rarity condition $m(\infty)<u$ is satisfied for all $u>0$ , as we have $0.9-1=-0.1<0<0.9\,u$ .

First we focus on the evaluation of the large-deviation asymptotics. For $s\rightarrow \infty$ , the Legendre transform I(s) is decreasing (increasing) if $u>\tfrac{1}{9}$ (if $u<\tfrac{1}{9}$ , respectively). For illustrational purposes we have plotted the functions $\alpha^\star(s)$ and I(s) in Figure 1, as a function of time s, for $u=5$ as well as $u=0.1$ . In the first instance, with $u=5$ , the function $I(\cdot)$ is decreasing, so the optimal $t^\star=\infty$ , whereas for $u=0.1$ we see that $I(\cdot)$ attains a minimal value at $t^\star=2.3$ .

Figure 1: The Legendre transform I(s) and the underlying optimal $\alpha^\star(s)$ parameter as a function of time s (for $s\in[0,5]$ ): (a,b) $u=5$ , (c,d) $u=0.1$ .

In Figure 2 we present, for different values of the initial number of obligors n and $u=5$ , the ruin probabilities as a function of time. This has been done relying on the iterative approach presented in Section 5.1. The double integral involved has been evaluated analytically for $n=1,2$ while numerical integration methods have been employed for $n>2$ . We observe that the ruin probability increases in the length of the time interval, as desired. The upper bound (as derived in Section 3.4) in this instance is given by 0.6065, and is independent of the number of obligors n. As can be observed, this upper bound is rather conservative, in particular when there are only a few obligors in the system.

Figure 2: Ruin probabilities over time: $p_n(u,t)$ as a function of t, for $n=1$ (bottom line) to $n=10$ (top line), with $u=5$ .

In a next experiment we study the performance of the importance sampling technique that was presented in Section 3.3. Figure 3(a) shows, for the initial capital reserve u being equal to 5, the estimates of the ruin probability as a function of time, obtained by simulation, using our importance sampling algorithm. The values nearly coincide with what is obtained by applying the naïve, direct simulation approach (i.e. without a change of measure); from Figure 2 we also observe that there is a highly accurate match with the values computed using the iterative approach of Section 5.1. Regarding the importance sampling simulations it is noted that we let the events $\mathscr{E}_j$ correspond to the event where the net cumulative loss process exceeds the initial level u (instead of nu), as u in this example corresponds to the unscaled initial capital level. The fact that we have used as many as $10^6$ runs guarantees estimates with a high precision. The importance-sampling-based approach substantially outperforms direct simulation, in that it greatly reduces the variance of the estimator, as can be observed in Figure 3(b,c).

Figure 3: (a) Ruin probabilities, as simulated by importance sampling: $p_n(u,t)$ as a function of time t. (b) Variance of the estimator under direct simulation as a function of t. (c) Variance of the estimator under importance sampling as a function of t. In all experiments we took $u=5$ .