2. Lévy-driven contagion models

In this section we construct a new framework for modelling event arrivals with contagion effects based on Lévy processes. That is, the intensity of the point process is set up to be a non-Gaussian OU process driven by a Lévy subordinator in cooperation with extra self-exciting jumps. Let us first define a simpler version without the self-exciting component.

Definition 2.1.

(Jump process with non-Gaussian intensity) A jump process with non-Gaussian intensity is a point process $N\equiv \{T_i\}_{i=1,2,\ldots}$ , i.e. $N_t = \sum_{i\geq 1}\textbf{1}_{\{ T_i \leq t\}}$ with the stochastic intensity $\lambda_t$ satisfying the stochastic differential equation

(2.1) $$\begin{equation} \textrm{d}\lambda_t = -\delta \lambda_t \,\textrm{d}t + \varrho\,\textrm{d}Z_t, \quad t\geq 0, \label{eqn1} \end{equation}$$

where

• $\varrho>0$ is an arbitrary positive constant,

• $\delta>0$ is the constant rate of exponential decay,

• $Z_t$ , with $Z_0=0$ , is a Lévy subordinator, which is called the background driving Lévy process (BDLP) of a non-Gaussian OU process.

This is a special case of Cox point processes. A slightly mathematical generalisation, though very useful for applications, is to further incorporate a feedback mechanism in the framework by adding a series of self-exciting jumps, i.e. simultaneous jumps (or ‘co-jumps’), in the point process and its intensity process. More precisely, this new framework, as a generalised version of the jump process of Definition 2.1, is defined via the stochastic intensity representation, as follows.

Definition 2.2.

(Self-exciting jump process with non-Gaussian intensity) $N_t$ is a self-exciting jump process with non-Gaussian intensity, if the intensity process of (2.1) is replaced by

(2.2) $$\begin{equation} \textrm{d}\lambda_t = -\delta \lambda_t \,\textrm{d}t + \varrho\,\textrm{d}Z_{t} + \textrm{d}J_t, \quad t\geq 0, \label{eqn2} \end{equation}$$

where the extra component $J_t$ is a pure-jump process specified by

(2.3) $$\begin{equation} J_t\,{:\!=}\, \sum^{N_t}_{i=1}X_i, \label{eqn3} \end{equation}$$

and $\{X_i\}_{i=1,2,\ldots}$ are the sizes of self-exciting jumps, as the expression (2.3) reveals that the jumps simultaneously occur in the point process $N_t$ and its intensity $\lambda_t$ , and hence the arrivals of jumps trigger more jumps afterwards, with $\mathcal{F}_{T^{-}_i}$ -measurable distribution function G(z), $z>0$ , occurring at the associated (ordered) arrival times $\{T_i\}_{i=1,2,\ldots}$ , respectively, and they are independent of $Z_t$ . It means that the functional form of the distribution function G(z) is revealed just before the arrival time $T_i$ . This distribution could have a highly general dependence structure $G(z)=G(z\!{\mid}\!{\cdot})$ . For example, it could depend on the initial intensity $\lambda_0$ , the past history of intensity at or just before the jump arrival times $\{T_k\}_{k=1,2,\ldots,i}$ , all past jump sizes $ \{X_k\}_{k=1,2,\ldots,i-1}$ , or the cumulated number of jumps $N_t$ , and so on, as long as we can record this information, for example,

$$G(z) = G ( z {\!{\mid}\!} T_1, T_2,\ldots,T_i, \lambda_0, \lambda_{T_1^{-}},\ldots, \lambda_{T_i^{-}}, \lambda_{T_1},\ldots, \lambda_{T_{i-1}} ).$$

This is similar to the adaptive model setting of [Reference Giesecke, Kakavand and Mousavi63], but it is different from the classical Hawkes process.

Similarly to the Hawkes process [Reference Hawkes and Oakes72], $N_t$ in Definition 2.2 can be equivalently redefined as a branching process via a cluster process presentation [Reference Daley and Vere-Jones43]. More precisely, $N_t$ is a cluster point process which consists of two types of points: immigrants and their offspring. The arrivals of immigrants follow a Cox process with non-Gaussian OU intensity (2.1). Each immigrant generates its offspring, each offspring would generate further offspring, and so on. The generation of any offspring follows a Cox process with exponentially decaying intensity $X^* \,{{\rm e}}^{-\delta (t-T^* )}$ , where $X^* {\,{\buildrel \mathcal{D} \over =}\,} X_i$ and T ^* is the arrival (birth) time of its previous generation. The superposition in [Reference Daley and Vere-Jones43] of all of these points forms our new self-exciting point process $N_t$ with the stochastic intensity (2.2).

Note that, given the initial intensity level $\lambda_0>0$ , the intensity process (2.2) can be equivalently expressed as

$$ \lambda_t = \underset{\text{{exogenous commonly shared} risk}}{\underbrace{ \lambda_0 \,{{\rm e}}^{-\delta t} + \varrho \int_0^t \,{{\rm e}}^{-\delta (t-s)} \,\textrm{d}Z_s }} + \underset{\text{{endogenous contagion risk}}}{\underbrace{\sum_{i=1}^{N_t}\,{{\rm e}}^{-\delta(t-T_i)} X_i, }} \quad t\geq 0, $$

which is positive and càdlàg. In fact, this new framework integrates two major types of risk sources. For example, in the context of credit risk or systemic risk, the first part (i.e. the first two terms) is to model the cyclical dependence of companies on some exogenous risk (e.g. movements of interest or FX rates) commonly shared in the entire market, and the cyclical oscillation is captured by the mean-reverting non-Gaussian OU process; fundamental common shocks are captured by the pure-jump process $Z_t$ . The second part (i.e. the last term) is to model the endogenous contagion risk due to the local interaction of companies in their business network, without which the overall risk would be underestimated.

The interarrival intensity process $\{\lambda_t \}_{T_i\le t < T_{i+1}}$ , for modelling exogenous commonly shared risk, is defined as the parts of intensity process excluding self-exciting jumps, i.e. (2.1), or

$$ \lambda_t = {{\rm e}}^{-\delta t} \lambda_0 + \varrho \int _0 ^ t {{\rm e}}^ {-\delta (t-s) } \,\textrm{d}Z_s, \quad t \in [T_i, T_{i+1}). $$

For instance, a sample path of the interarrival intensity process (without self-exciting component) within the time period $t \in[0, 5]$ when the BDLP is a gamma process is represented in Figure 1.

Fig. 1.

A sample path of interarrival intensity process (without self-exciting component) within the time period $t \in[0, 5]$ when the BDLP is a gamma process.

This framework is the generalisation of several classical models in the literature. If there are no self-exciting jumps, i.e. $X_i\equiv0$ for any i, then the intensity process (2.2) is a non-Gaussian OU process [Reference Barndorff-Nielsen, Shephard, Barndorff-Nielsen, Resnick and Mikosch14, Reference Barndorff-Nielsen and Shephard15, Reference Barndorff-Nielsen and Shephard16, Reference Barndorff-Nielsen and Shephard17, Reference Barndorff-Nielsen, Jensen and Sørensen19]. If there is no BDLP $Z_t$ , then the point process $N_t$ is a generalised Hawkes process [Reference Hawkes70, Reference Hawkes71] with random marks. If the BDLP $Z_t$ is trivially a subordinator of compound Poisson, then $N_t$ is a dynamic contagion process [Reference Dassios and Zhao46].

For notational simplicity, we denote the Lévy measure of BDLP $Z_t$ by $\nu$ , the associated Laplace exponent and mean at unit time by

$$ \Phi(u)\,{:\!=}\, \int_0^\infty ( 1 - {{\rm e}}^{-u y } ) \nu({\rm d}y), \quad \mu_Z \,{:\!=}\, \mathbb{E}[Z_1] = \int_0^\infty y\nu({\rm d}y), \quad u>0, $$

the Laplace transform, mean of self-exciting jump sizes and a constant, respectively, by

$$ \hat{g}(u) \,{:\!=}\, \int_0^\infty {{\rm e}}^{-uy} \,\textrm{d}G(y), \quad \mu_{G} \,{:\!=}\, \int_0^\infty y \,\textrm{d}G(y), \quad \eta\,{:\!=}\,\delta-\mu_{G}. $$

They are all assumed to be finite. In addition, we denote the $(i+1)$ th interarrival time by

$$ \tau _ {i+1} \,{:\!=}\, T _ {i + 1} - T_i, \quad i = 0,1,2,\ldots, \quad T_0 = 0, $$

and the cumulative intensity process at time t by $\Lambda_t \,{:\!=}\, \int _ 0 ^ t \lambda_u \,{\rm d}u$ .

One may be interested in their basic distributional properties such as means and Laplace transforms, and we provide some brief results for them in Appendix A with proofs. In particular, the conditional expectation of point process is provided here in Proposition 2.1, as it will be used later as a simple and general benchmark for numerically validating our newly developed simulation algorithms.

Proposition 2.1.

(Conditional expectation of point process.) The expectation of $N_{t+s}$ conditional on $N_t$ and $\lambda_t$ is given by

$$ \begin{equation}\mathbb{E} [ N_{t+s} {\!{\mid}\!} N_t, \lambda_t ] = \begin{cases} N_t + \dfrac{\varrho \mu_Z }{\eta} s + \bigg( \lambda_t -\dfrac{\varrho \mu_Z}{\eta} \bigg) \dfrac{1- {{\rm e}}^{-\eta s}}{\eta} & \eta\ne 0,\\[14pt] N_t + \lambda_t s + \dfrac{1}{2}\varrho\mu_Z s^2 & \eta= 0, \end{cases} \quad s>0.\end{equation}$$

3. General framework for exact simulation

In this section we outline an exact simulation framework based on exact distributional decomposition for a general point process of Lévy-driven non-Gaussian OU intensity with and without self-exciting jumps as defined in Definition 2.2 and 2.1, respectively. The entire simulation scheme can be decomposed into three major steps.

Conditional on the current arrival time $T_i$ and the associated intensity level $\lambda_{T_i}$ , generate the next interarrival time $\tau_{i+1}$ .

Further conditional on the realisation of this interarrival time $\tau_{i+1}$ , generate the pre-jump intensity level $\lambda_{T_i+\tau_{i+1}^-}$ right before the next arrival time $T_{i+1}=T_i+\tau_{i+1}$ .

Add a self-exciting jump size $X_{i+1}$ upon the intensity process $\lambda_t$ and one unit in the point process $N_t$ both at the next arrival time $T_{i+1}=T_i+\tau_{i+1}$ .

By recursively implementing the three steps above, a full path of the point process $N_t$ in any time horizon can be exactly produced without bias. A graphical illustration of this proposed algorithm design is provided in Figure 2.

Fig. 2.

Exact simulation procedures for a path of point process $N_t$ and the skeleton of its intensity process $\lambda_t$ around the period $[T_i, T_{i+1}]$ .

The third step is in fact straightforward. In particular, if $X_{i+1} \equiv 0$ for any i, then it corresponds to the version without self-exciting jumps. In order to execute the first two steps, we have first to further investigate the joint distributional properties of the next interarrival time $\tau_{i+1}$ and the next pre-jump intensity level $\lambda_{T_i+\tau_{i+1}^-}$ , which can be characterised by the conditional joint transform as follows.

Theorem 3.1.

(Joint transform of pre-jump intensity and cumulative intensity.) Conditional on the intensity level $\lambda_{T_i}$ at the ith arrival time $T_i$ , the joint transform of $ (\lambda_{T_i+\tau^-}, \Lambda_{T_i+\tau}-\Lambda_{T_i})$ for any given period $\tau \in (0, \tau_{i+1} )$ is given by

(3.1) $$\begin{align} & \mathbb{E} [ {{\rm e}}^{-v\lambda_{T_i+\tau^-}} \,{{\rm e}}^{- (\Lambda_{T_i+ \tau }-\Lambda_{T_i})} {\!{\mid}\!} \lambda_{T_i} ] \notag \\ &\qquad = \exp \bigg( -\bigg[ \dfrac{1}{\delta}+ \bigg( v - \dfrac{1}{\delta}\bigg)w\bigg] \lambda_{T_i} - \varrho\int_v^{\frac{1}{\delta}+(v-\frac{1}{\delta} )w} \dfrac{\Phi(u)}{1-\delta u} \,{\rm d}u\bigg), \quad \tau \in (0, \tau_{i+1} ), \label{eqn4}\end{align}$$

where $w \,{:\!=}\, {{\rm e}}^{ -\delta \tau }$ .

Proof. Given the ith arrival time $T_i$ , the infinitesimal generator of ${(\Lambda,\lambda, t)}$ within the period $t\in [T_i, T_i+\tau_{i+1})$ acting on any function $f(\Lambda,\lambda,t)$ within its domain $\Omega(\mathcal{A})$ is given by

(3.2) $$\begin{equation} \mathcal{A}f(\Lambda,\lambda,t) = \dfrac{\partial f}{\partial t} - \delta\lambda\dfrac{\partial f}{\partial \lambda}+\lambda\dfrac{\partial f}{\partial \Lambda} + \varrho\bigg\{ \int_0^\infty [f(\Lambda,\lambda+y,t) - f(\Lambda,\lambda,t)]\nu({\rm d}y) \bigg\}. \label{eqn5} \end{equation}$$

Consider a function

(3.3) $$\begin{equation} f(\Lambda,\lambda,t) = {{\rm e}}^{-\tilde{v} \Lambda} \,{{\rm e}}^{-\lambda A(t)} \,{{\rm e}}^{R(t)}, \quad \tilde{v} \in \mathbb{R}^+, \label{eqn6} \end{equation}$$

where A(t) and R(t) are deterministic and differentiable functions with respect to t. Substituting (3.3) into (3.2) and setting $\mathcal{A}f=0$ , we have

$$ -\lambda A'(t) +R'(t) +\delta \lambda A(t) -\tilde{v}\lambda- \varrho\Phi (A(t)) = 0. $$

Since this equation holds for any $\lambda$ and $\Lambda$ , it is equivalent to the equations

$$ A'(t) = \delta A(t) -\tilde{v}, \quad R'(t) = \varrho\Phi (A(t)). $$

Hence, for any time $t\in [T_i, T_i+\tau_{i+1} )$ , we have

$$ A(t) = k \,{{\rm e}}^{ \delta t} -\tilde{v} \dfrac{ {{\rm e}}^{\delta t} -1 }{\delta}, \quad R(t) = \varrho\int_0^t \Phi\bigg( k \,{{\rm e}}^{\delta s} - \tilde{v} \dfrac{{{\rm e}}^{\delta s} -1 }{\delta} \bigg) \,{\rm d}s, \quad k\in\mathbb{R}^+. $$

By the basic property of infinitesimal generator in [Reference Dassios and Embrechts45], we have the martingale

$$\begin{equation}{{\rm e}}^{-\tilde{ v} \Lambda_t} \exp\bigg( -\bigg( k \,{{\rm e}}^{\delta t} - \tilde{v} \dfrac{{{\rm e}}^{\delta t} -1 }{\delta} \bigg) \lambda_t + \varrho\int_0^t \Phi\bigg( k \,{{\rm e}}^{\delta s} - \tilde{v} \dfrac{{{\rm e}}^{\delta s } -1}{\delta} \bigg) \,{\rm d}s \bigg).\end{equation}$$

Setting $\tilde{v} = 1$ and $A (T_i+\tau^- )= v$ for any $\tau \in (0, \tau_{i+1} )$ and using the martingale property, we have

$$\begin{align*} & \mathbb{E} [ {{\rm e}}^ {-v\lambda_{T_i+\tau^-}} \,{{\rm e}}^ {- (\Lambda_{T_i+ \tau}-\Lambda_{T_i} )} {\!{\mid}\!} \lambda_{T_i}] \\ &\qquad = \exp \bigg( -\bigg[ \dfrac{1}{\delta}+ \bigg(v - \dfrac{1}{\delta} \bigg)w \bigg] \lambda_{T_i} - \varrho\int _ 0 ^ {\tau} \Phi\bigg( \dfrac{1}{\delta}+ \bigg( v - \dfrac{ 1}{\delta} \bigg) \,{{\rm e}}^ {-\delta s}\bigg) \,{\rm d}s\bigg). \end{align*}$$

Hence, we can immediately obtain (3.1) by the change of variable $u = \frac{1}{\delta}+ ( v - \frac{1}{\delta}) \,{{\rm e}}^ {-\delta s}$ .

Theorem 3.1 provides us with a crucial tool for further investigating the distributional properties of the interarrival time $\tau_{i+1}$ and the pre-jump intensity level $\lambda_{T_i+\tau_{i+1}^-}$ , jointly and separately, which later leads to their efficient algorithms for exact simulation as follows. Note that, however, the main mathematical contribution of this paper is the exact distributional decomposition for our new point processes via Theorem 3.1 rather than deriving the transform in Theorem 3.1.

3.1. Exact simulation of interarrival time

Let us first outline how to simulate the interarrival time. Given the intensity level $\lambda_{T_i}$ at the ith arrival time $T_i$ , interestingly, the $(i+1)$ th interarrival time $\tau_{i+1}$ can be exactly expressed as the minimum of two much simpler random variables $V^*_{T_i}$ and V ^* , where

$V^*_{T_i}$ is a defective random variable, which can be directly generated by an explicit inverse transform,

V ^* is a well-defined random variable, which can be exactly simulated by a simplified version of the classical thinning scheme [Reference Lewis and Shedler79].

Algorithm 3.1.

(Exact simulation of interarrival time.) Conditional on the intensity level $\lambda_{T_i}$ , the next interarrival time $\tau_{i+1}$ can be exactly simulated via

$$\begin{equation}\tau_{i+1}{\,{\buildrel \mathcal{D} \over =}\,} V^* \wedge V^*_{T_i} D_i > 0, \\ {V^* D_i < 0},\end{equation}$$

where

• $D_i$ is simulated via

  $$ D_i: = 1 + \dfrac{\delta}{\lambda_{T_i}}\ln U_1, \quad U_1 \sim \operatorname{\textsf{U}} [0,1]{,} $$

• $V^*_{T_i}$ is a simple defective random variable with

  $ \mathbb{P } \{ V^*_{T_i} =\infty \} = \exp ( -\tfrac{1}{\delta} \lambda_{T_i} )$ ,

  (3.4) $$\begin{equation} V^*_{T_i} {\,{\buildrel \mathcal{D} \over =}\,} - \dfrac{1}{\delta} \ln D_i, \quad D_i>0 {,} \label{eqn7} \end{equation}$$

• V ^* is the first arrival time of a non-homogeneous Poisson process with the rate function

  (3.5) $$\begin{equation} \zeta_t\,{:\!=}\, \varrho \Phi ( G_0(t)), \quad G_0(u)\,{:\!=}\, \dfrac{1- {{\rm e}}^ {-\delta u}}{\delta}, \quad u \ge 0, \label{eqn8} \end{equation}$$
  and it can be exactly simulated via the simplified thinning scheme of Algorithm 3.2.

Proof.

Setting $v=0$ in (3.1) of Theorem 2, we have

$$ \mathbb{P } \{ \tau_{i+1} > \tau \!{\mid}\! \lambda_{T_i} \} = \mathbb{E} [ {{\rm e}}^ {- ( \Lambda_{T_i+\tau} - \Lambda_{T_i} )} \!{\mid}\! \lambda_{T_i}] = \mathbb{P } \{ V^*> \tau \} \times \mathbb{P } \{ V_{T_i}^* > \tau \}, $$

where

(3.6) $$\begin{equation} \mathbb{P } \{ V^*> \tau \} = \exp \bigg( - \varrho \int _ 0 ^ \tau \Phi ( G_0(u)) \,{\rm d}u \bigg), \quad \mathbb{P } \{ V^*_{T_i}> \tau \} = {{\rm e}}^{ - G_0(\tau) \lambda_{T_i}}. \label{eqn9} \end{equation}$$

This implies that the next interarrival time $\tau_{i+1}$ conditional on the current intensity level $\lambda_{T_i}$ can be expressed as the minimum of two independent random variables V ^* and $ V^*_{T_i}$ . Note that $V^*_{T_i}$ is a defective random variable that has a mass probability at the infinity, since the CDF of $V^*_{T_i}$ is

$$ F_{V^*_{T_i}}(\tau) = 1 - {{\rm e}}^{ - G_0(\tau) \lambda_{T_i} }, $$

with

$$ F_{V^*_{T_i}}(\infty)=1- \exp\bigg(-\dfrac{1}{\delta}\lambda_{T_i}\bigg)<1, $$

and the density $f_{V^*_{T_i}}(\tau)>0$ for any $\tau>0$ . Obviously, if $D_i>0$ , then $V^*_{T_i}$ can be exactly simulated using the explicit inverse transform (3.4), whereas V ^* can be interpreted as the first arrival time from a non-homogeneous Poisson process, and it can be exactly simulated via Algorithm 3.2 as follows.

Algorithm 3.2

(Simplified thinning scheme.) V ^* can be exactly simulated by the following steps.

1. 1. Initialise the candidate time $\skew3\tilde{t}=0$ .

2. 2. Generate an exponential random variable $E^* \sim \operatorname{\textsf{Exp}}( \zeta_\infty )$ where

   (3.7) $$\begin{equation} \zeta_\infty \,{:\!=}\, \lim_{t \rightarrow \infty} \zeta_t = \varrho \Phi\bigg( \dfrac{1}{\delta} \bigg), \label{eqn10} \end{equation}$$
   and set $\skew3\tilde{t} = \skew3\tilde{t} + E^*$ .

3. 3. Generate a uniform random variable $U_2 \sim \operatorname{\textsf{U}} [0,1]$ ,

   • if $U_2 \leq \zeta_{\skew3\tilde{t}}/\zeta_\infty$ , then accept this candidate by setting $V^* = \skew3\tilde{t}$ ,

   • if $U_2 > \zeta_{\skew3\tilde{t}}/\zeta_\infty$ , then reject this candidate, and go back to step 2 and continue.

Proof.

Since $\zeta_t$ in (3.5) is a strictly increasing and concave function of time t with the initial value $\zeta _0 = 0$ at time $t=0$ , the maximum level is $\zeta_\infty$ . Then the algorithm above is in fact a simplified version of the classical thinning scheme [Reference Lewis and Shedler79] where only the first arrival time within the period of [0, t] is recorded.

3.2. Exact simulation of pre-jump intensity level

Conditional on the realisation of interarrival time $\tau_{i+1}$ as generated by Algorithm 3.1, the Laplace transform of the next pre-jump intensity level $\lambda_{T_i+\tau_{i+1}^-}$ is provided as follows.

Theorem 3.2.

(Laplace transform of pre-jump intensity.) Conditional on the intensity level $\lambda_{T_i}$ and the $(i+1)$ th interarrival time $\tau_{i+1}$ , the Laplace transform of pre-jump intensity level $\lambda_{T_i+\tau_{i+1}^-}$ is given by

(3.8) $$\begin{align} &\mathbb{E} [ {{\rm e}}^ {-v \lambda_{T_i+\tau_{i+1}^-}} \!{\mid}\! \tau_{i+1}=\tau, \lambda_{T_i}] \notag \\[6pt] &\qquad = {{\rm e}}^{{-vw\lambda_{T_i}}} \times \exp\bigg(-\varrho \int_{v}^{\frac{1}{\delta}+ (v-\frac{1}{\delta})w}\dfrac{\Phi(u)-\Phi\big( u-\tfrac{u-\frac{1}{\delta}}{v-\frac{1}{\delta}}v \big)}{1-\delta u}\,{\rm d}u\bigg) \notag \\[6pt] &\qquad \quad\, \times \dfrac{\displaystyle \dfrac{\varrho}{\delta} \int_{0}^{\infty}\,{{\rm e}}^{-vs}\int_{s}^{\frac{s}{w}}\,{{\rm e}}^{\frac{s}{\delta}}\,{{\rm e}}^{-\frac{y}{\delta}}\nu({\rm d}y)\,{\rm d}s+w\lambda_{T_i}} { \displaystyle \dfrac{\varrho}{\delta} \int_{0}^{\infty}\int_{s}^{\frac{s}{w}} \,{{\rm e}}^{\frac{s}{\delta}}\,{{\rm e}}^{-\frac{y}{\delta}}\nu({\rm d}y)\,{\rm d}s+w\lambda_{T_i}}. \label{eqn11}\end{align}$$

Proof. Note that the density function of the $(i+1)$ th interarrival time conditional on the intensity level $\lambda_{T_i}$ is

$$ \mathbb{P } \{ \tau_{i+1} \in \textrm{d}\tau \!{\mid}\! \lambda_{T_i} \} = \mathbb{E}\bigg[ \lambda_{T_i+\tau^-}\exp\bigg( -\int _ {T_i} ^ {T_i+\tau}\lambda_u \,{\rm d}u \bigg) \Bigm| \lambda_{T_i}\bigg]\,\textrm{d}\tau, $$

which implies that

$$ \mathbb{E}\big[ {{\rm e}}^{-v\lambda_{T_i+\tau_{i+1}^-}}\textbf{1}_{ \{\tau_{i+1}\in \textrm{d}\tau \}} \!{\mid}\! \lambda_{T_i} \big] = \mathbb{E}\bigg[ \lambda_{T_i+\tau^-}\,{{\rm e}}^{-v\lambda_{T_i+\tau^-}} \exp\bigg(- \int_{T_i} ^ {T_i+\tau}\lambda_u \,{\rm d}u \bigg) \Bigm| \lambda_{T_i}\bigg]\,\textrm{d}\tau. $$

Hence, we have

(3.9) $$\begin{equation} \mathbb{E}\big[ {{\rm e}}^{-v\lambda_{T_i+\tau_{i+1}^-}} \!{\mid}\! \tau_{i+1}=\tau, \lambda_{T_i}\big] = \dfrac{\mathbb{E}\big[\lambda_{T_i+\tau^-}\,{{\rm e}}^{-v\lambda_{T_i+\tau^-}}\,{{\rm e}}^{- (\Lambda_{T+\tau}-\Lambda_{T_i} )} \!{\mid}\! \lambda_{T_i}\big]} {\mathbb{E}\big[\lambda_{T_i+\tau^-}\,{{\rm e}}^{- (\Lambda_{T_i+\tau}-\Lambda_{T_i})} \!{\mid}\! \lambda_{T_i}\big]} . \label{eqn12} \end{equation}$$

The numerator of (3.9) can be obtained by differentiating the joint transform (3.1) with respect to v, that is,

$$\begin{align*} & \mathbb{E}\big[\lambda_{T_i+\tau^-}\,{{\rm e}}^{-v\lambda_{T_i+\tau^-}}\,{{\rm e}}^{- (\Lambda_{T_i+\tau}-\Lambda_{T_i} )} \!{\mid}\! \lambda_{T_i}\big] \\[6pt] & \qquad = -\dfrac{\partial }{\partial v} \mathbb{E} \big[ {{\rm e}}^{- v \lambda_{T_i+\tau^-}} \,{{\rm e}}^{- (\Lambda_{T_i+\tau^-}-\Lambda_{T_i} )} \!{\mid}\! \lambda_{T_i} \big] \notag \\[6pt] & \qquad = -\dfrac{\partial}{\partial v}\bigg[\exp\bigg(-\bigg[\dfrac{1}{\delta}+ \bigg(v-\dfrac{1}{\delta} \bigg)\,{{\rm e}}^{-\delta \tau }\bigg]\lambda_{T_i}\bigg) \\[6pt] &\qquad\quad\, \times \exp\bigg(-\varrho\int_{0}^{\tau}\Phi\bigg(\dfrac{1}{\delta}+\bigg(v-\dfrac{1}{\delta}\bigg)\,{{\rm e}}^{-\delta u}\bigg)\,{\rm d}u\bigg)\bigg] \notag \\[6pt] & \qquad = -\dfrac{\partial}{\partial v}\big[ {{\rm e}}^{{-G_v(\tau)\lambda_{T_i}}} \times {{\rm e}}^{{-\varrho F_v(\tau)}} \big] \notag \\[6pt] & \qquad = -\bigg[\varrho\dfrac{\partial}{\partial v}F_v(\tau) + \dfrac{\partial}{\partial v} G_v(\tau) \bigg] \times {{\rm e}}^{-G_v(\tau)\lambda_{T_i}} \,{{\rm e}}^{-\varrho F_v(\tau)}, \notag \end{align*}$$

where

$$ G_v(u) \,{:\!=}\, \dfrac{1}{\delta}+\bigg(v-\dfrac{1}{\delta}\bigg)\,{{\rm e}}^{-\delta u }, \quad F_v(\tau) \,{:\!=}\, \int_{0}^{\tau}\Phi (G_v(u))\,{\rm d}u, $$

and

$$\begin{align*} \dfrac{\partial }{\partial v} G_v(\tau) &= -\lambda_{T_i} \,{{\rm e}}^{-\delta \tau}, \\[5pt] \dfrac{\partial}{\partial v}F_v(\tau) &= \int_{0}^{\tau}\int_{0}^{\infty} y\,{{\rm e}}^{-\delta u}\,{{\rm e}}^{- [\frac{1}{\delta}+ (v-\frac{1}{\delta})\,{{\rm e}}^{-\delta u} ]y}\nu({\rm d}y)\,{\rm d}u. \end{align*}$$

Note that $\nu$ is the Lévy measure for a general BDLP $Z_t$ , so we have

(3.10) $$\begin{align} & \mathbb{E}\big[ {{\rm e}}^{-v\lambda_{T_i+\tau_{i+1}^-}} \!{\mid}\! \tau_{i+1}=\tau, \lambda_{T_i}\big] \notag \\[5pt] &\qquad = \dfrac{\displaystyle\bigg(\varrho\int_{0}^{\tau}\int_{0}^{\infty}y\,{{\rm e}}^{-\delta u}\,{{\rm e}}^{- [\frac{1}{\delta}+ (v-\frac{1}{\delta})\,{{\rm e}}^{-\delta u} ]y}\nu({\rm d}y)\,{\rm d}u+w\lambda_{T_i}\bigg) \times {{\rm e}}^{-G_v(\tau)\lambda_{T_i}} \,{{\rm e}}^{-\varrho F_v(\tau)}}{\displaystyle\bigg(\varrho\int_{0}^{\tau}\int_{0}^{\infty}y\,{{\rm e}}^{-\delta u}\,{{\rm e}}^{-\frac{1}{\delta} (1-{{\rm e}}^{-\delta u})y}\nu({\rm d}y)\,{\rm d}u+w\lambda_{T_i}\bigg)\times {{\rm e}}^{-G_0(\tau)\lambda_{T_i}} \,{{\rm e}}^{-\varrho F_0(\tau)}} \notag \\[5pt] &\qquad = \dfrac{\displaystyle\varrho\int_{0}^{\tau}\int_{0}^{\infty}y\,{{\rm e}}^{-\delta u}\,{{\rm e}}^{- [\frac{1}{\delta}+ (v-\frac{1}{\delta})\,{{\rm e}}^{-\delta u} ]y}\nu({\rm d}y)\,{\rm d}u+w \lambda_{T_i}}{\displaystyle\varrho\int_{0}^{\tau}\int_{0}^{\infty}y\,{{\rm e}}^{-\delta u}\,{{\rm e}}^{-\frac{1}{\delta} (1-{{\rm e}}^{-\delta u})y}\nu({\rm d}y)\,{\rm d}u+w \lambda_{T_i}} \notag \\[5pt] &\qquad\quad\, \times \dfrac{\mathbb{E}\big[{{\rm e}}^{-v\lambda_{T_i+\tau^-}}\,{{\rm e}}^{- (\Lambda_{T_i+\tau}-\Lambda_{T_i})} \!{\mid}\! \lambda_{T_i}\big]}{\mathbb{E}\big[{{\rm e}}^{- (\Lambda_{T_i+\tau}-\Lambda_{T_i})} \!{\mid}\! \lambda_{T_i}\big]}. \label{eqn13}\end{align}$$

The first term of (3.10) can be calculated more explicitly as

(3.11) $$\begin{align} & \dfrac{\displaystyle\varrho\int_{0}^{\tau}\int_{0}^{\infty}y\,{{\rm e}}^{-\delta u}\,{{\rm e}}^{- [\frac{1}{\delta}+ (v-\frac{1}{\delta})\,{{\rm e}}^{-\delta u}]y}\nu({\rm d}y)\,{\rm d}u+w\lambda_{T_i}} {\displaystyle\varrho\int_{0}^{\tau}\int_{0}^{\infty}y\,{{\rm e}}^{-\delta u}\,{{\rm e}}^{-\frac{1}{\delta} (1-{{\rm e}}^{-\delta u})y}\nu({\rm d}y)\,{\rm d}u+w\lambda_{T_i}}\notag \\[5pt] &\qquad = \dfrac{\displaystyle\varrho\int_{v}^{\frac{1}{\delta}+ (v-\frac{1}{\delta})w} \int_{0}^{\infty}y\dfrac{u-\frac{1}{\delta}}{v-\frac{1}{\delta}}\,{{\rm e}}^{-uy}\nu({\rm d}y)\dfrac{{\rm d}u}{1-\delta u}+w\lambda_{T_i}}{\displaystyle\varrho\int_{v}^{\frac{1}{\delta}+ (v-\frac{1}{\delta})w} \int_{0}^{\infty}y\dfrac{u-\frac{1}{\delta}}{v-\frac{1}{\delta}}\,{{\rm e}}^{- (u-\frac{u-\frac{1}{\delta}}{v-\frac{1}{\delta}}v )y}\nu({\rm d}y)\dfrac{{\rm d}u}{1-\delta u}+w\lambda_{T_i}} \notag \\[5pt] & \qquad = \dfrac{ \displaystyle\dfrac{\varrho}{1-\delta v} \int_{v}^{\frac{1}{\delta}+ (v-\frac{1}{\delta} )w}\int_{0}^{\infty}y\,{{\rm e}}^{- (u-\frac{1}{\delta})y}\,{{\rm e}}^{-\frac{y}{\delta}}\nu({\rm d}y)\,{\rm d}u+w\lambda_{T_i}} { \displaystyle\dfrac{\varrho}{1-\delta v} \int_{v}^{\frac{1}{\delta}+ (v-\frac{1}{\delta})w}\int_{0}^{\infty}y\,{{\rm e}}^{-(u-\big(\frac{1}{\delta}-\tfrac{u-\frac{1}{\delta}}{v-\frac{1}{\delta}}v\big)y}\,{{\rm e}}^{-\frac{y}{\delta}}\nu({\rm d}y)\,{\rm d}u+w\lambda_{T_i}}\notag \\[5pt] & \qquad = \dfrac{\displaystyle -\dfrac{\varrho}{\delta} \int_{1-\delta v}^{(1-\delta v)w}\int_{0}^{\infty}y\,{{\rm e}}^{\frac{((1-\delta v)zy)}{\delta}}\,{{\rm e}}^{-\frac{y}{\delta}}\nu({\rm d}y)\,\textrm{d}z+w\lambda_{T_i}}{\displaystyle - \dfrac{\varrho}{\delta} \int_{1-\delta v}^{(1-\delta v)w}\int_{0}^{\infty}y\,{{\rm e}}^{\frac{zy}{\delta}}\,{{\rm e}}^{-\frac{y}{\delta}}\nu({\rm d}y)\,\textrm{d}z+w\lambda_{T_i}} \notag \\[5pt] & \qquad = \dfrac{ \displaystyle\dfrac{\varrho}{\delta} \int_{0}^{\infty}\,{{\rm e}}^{-vs}\int_{s}^{\frac{s}{w}} \,{{\rm e}}^{\frac{s}{\delta}}\,{{\rm e}}^{-\frac{y}{\delta}}\nu({\rm d}y)\,{\rm d}s+w\lambda_{T_i}} { \displaystyle\dfrac{\varrho}{\delta} \int_{0}^{\infty}\int_{s}^{\frac{s}{w}} \,{{\rm e}}^{\frac{s}{\delta}}\,{{\rm e}}^{-\frac{y}{\delta}}\nu({\rm d}y)\,{\rm d}s+w\lambda_{T_i}}. \label{eqn14}\end{align}$$

The denominator of the second term of (3.10) can also be obtained nicely by setting $v = 0$ in the joint transform (3.1), that is,

$$ \mathbb{E}\big[ {{\rm e}}^ {- (\Lambda_{T_i+\tau}-\Lambda_{T_i} )} \!{\mid}\! \lambda_{T_i} \big] = \exp \bigg( - \dfrac{1-w}{\delta} \lambda_{T_i} - \varrho \int_v^{\frac{1}{\delta}+ (v-\frac{1}{\delta})w} \dfrac{\Phi\big( u - \tfrac{u-\frac{1}{\delta}}{v-\frac{1}{\delta}} v\big)}{1-\delta u} \,{\rm d}s\bigg). $$

Therefore the second term of (3.10) can be expressed by

(3.12) $$\begin{align} & \dfrac{\mathbb{E}\big[ {{\rm e}}^{-v\lambda_{T_i+\tau^-}} \,{{\rm e}}^{- (\Lambda_{T_i+\tau}-\Lambda_{T_i} )} \!{\mid}\! \lambda_{T_i}\big] }{\mathbb{E}\big[{{\rm e}}^{- (\Lambda_{T_i+\tau}-\Lambda_{T_i})} \!{\mid}\! \lambda_{T_i} \big]} \notag \\[5pt] & \qquad = \dfrac{\exp \bigg( -\bigg[ \dfrac{1}{\delta}+ \bigg( v - \dfrac{1}{\delta} \bigg)w\bigg] \lambda_{T_i} \bigg) \exp\bigg(\displaystyle - \varrho \int_v ^ {\frac{1}{\delta}+ (v-\frac{1}{\delta})w} \dfrac{\Phi(u)}{1-\delta u} \,{\rm d}u\bigg)} {\exp \bigg( - \dfrac{1}{\delta}(1-w) \lambda_{T_i} \bigg) \exp\bigg(\displaystyle - \varrho\int_v ^{\frac{1}{\delta}+ (v-\frac{1}{\delta} )w} \dfrac{\Phi\big( u - \tfrac{u-\frac{1}{\delta}}{v-\frac{1}{\delta}} v\big)}{1-\delta u} \,{\rm d}u\bigg)} \notag \\[5pt] & \qquad = {{\rm e}}^{{-v w \lambda_{T_i}} } \times \exp\bigg( - \varrho \int_v ^{\frac{1}{\delta}+ (v-\frac{1}{\delta})w} \dfrac{\Phi(u) -\Phi\big(u - \tfrac{u-\frac{1}{\delta}}{v-\frac{1}{\delta}} v\big)}{1-\delta u} \,{\rm d}u\bigg). \label{eqn15}\end{align}$$

Finally, we obtain (3.8) immediately by combining the results from (3.11) and (3.12).

Apparently, given the ith arrival time $T_i$ and the $(i+1)$ th interarrival time $\tau_{i+1}$ , the pre-jump intensity level $\lambda_{T_i+\tau^-_{i+1}}$ can be simulated by the numerical inversion of the Laplace transform (3.8) for any Lévy-driven contagion process once the associated Lévy measure $\nu$ (and the Laplace exponent $\Phi$ ) are specified. Indeed, exact simulation for stochastic processes based on the numerical inversion of the Laplace or Fourier transforms has been widely adopted in the literature; see [Reference Broadie and Kaya24], Reference Glasserman and Liu66], [Reference Chen, Feng and Lin32], [Reference Cai, Song and Chen28], [Reference Kang, Kang and Lee74].

However, for some subclasses such as the very popular specifications of gamma and tempered stable BDLPs, quite remarkably, based on Theorem 3.2 the pre-jump intensity level can be exactly decomposed into several simple elements, each of which can easily be simulated exactly without any numerical inversion procedure. In fact, this exact decomposition approach appropriately breaks the Lévy measure of the subordinator, and thereby it can be achieved by developing an exact distributional decomposition through Laplace-transform representations. In this paper, our focus is mainly on this decomposition approach, as it leads to a very efficient simulation algorithm for exactly sampling the whole point process, and more importantly, it does not involve additional discretisation or truncation errors which are inevitable in the numerical inversion approach. We will present our discovery based on the decomposition approach in much more detail later in Section 4.

4. Typical examples: gamma and tempered stable contagion models

For model implementation, one needs to further specify the BDLP in an explicit form. Probably the most widely used and representative Lévy subordinators in the literature are the gamma process and the tempered stable (TS) subordinator. The term tempered stable in the context sometimes refers to exponential tilted stable (see [Reference Sato94] and [Reference Rosiski92] for more information). More precisely, they are two typical examples of Lévy processes with stationary, independent, and non-negative increments starting at 0 with finite variation and infinite activity, which means that there are infinitely many small jumps within any finite time interval (for some strong evidence from financial data see [Reference Li, Wells and Cindy80], [Reference At-Sahalia and Jacod3], [Reference At-Sahalia and Jacod4], Reference Lee and Hannig78]); these processes behave very differently from the trivial case of compound Poisson process. Because of this nature, they cannot be simulated exactly by traditional methods for discretising sample paths. We provide the definitions for them using the Lévy measures (4.1) and (4.3) in Definitions 4.1 and 4.2, respectively.

Definition 4.1.

(Gamma distribution and gamma process.) The gamma distribution with shape parameter a and rate parameter b, denoted by $\operatorname{\textsf{Gamma}}(a,b)$ , has the density function

$$f_{\operatorname{\textsf{Gamma}}(a,b)}(s) = \dfrac{b^a}{\Gamma(a)}s^{a-1} \,{{\rm e}}^{-b s}, \quad s> 0, \quad a,b>0,$$

where $\Gamma({\cdot})$ is the gamma function, i.e. $\Gamma(u)\,{:\!=}\, \int_0^\infty s^{u-1} \,{{\rm e}}^{-s}\,{\rm d}s$ . The gamma process $\{G_t, t\,{\ge}\, 0 \}$ is a pure-jump increasing Lévy process with independently gamma-distributed increments satisfying $G_1\sim \operatorname{\textsf{Gamma}}(a,b)$ , $G_0=0$ , and it has Lévy measure

(4.1) $$\begin{equation} \nu({\rm d}s) = a s^{-1} \,{{\rm e}}^{-b s} \,{\rm d}s, \label{eqn17} \end{equation}$$

and the Laplace exponent

(4.2) $$\begin{equation} \Phi(u)=a \ln \bigg( 1 + \dfrac{u}{b} \bigg), \label{eqn18} \end{equation}$$

and the mean at unit time

$$\begin{equation}\mu_Z=\mathbb{E}[Z_1]=\dfrac{a}{b}.\end{equation}$$

The gamma distribution, as a well-known right-skewed distribution, is approximated by a normal distribution when its mean is large, and it has many computational conveniences. A well-known convenience, similar to the normal distribution, is that if $\{Y_i\}_{i=1,2,\ldots,m}$ are independent, identically distributed (i.i.d.) random variables following a gamma distribution with mean $\mu_{Y}$ and variance $\sigma^2_{Y}$ , then the sum $\sum^m_{i=1}Y_i$ is gamma-distributed with mean $m \mu_{Y}$ and variance $m\sigma^2_Y$ . Thus it is a very useful building block as the risk driver. Therefore, we provide applications for the gamma-contagion model later in Section 6.

Definition 4.2.

(Tempered stable (TS) distribution and tempered stable (TS) subordinator.) The positive tempered stable (TS) distribution, abbreviated as $\operatorname{\textsf{TS}} (\alpha,\beta,\theta)$ , is defined by its Lévy measure

(4.3) $$\begin{equation}\nu({\rm d}y)=\dfrac{\theta}{y^{\alpha+1}}\,{{\rm e}}^{-\beta y}\,{\rm d}y,\quad y\geq 0,\quad \alpha\in(0,1), \quad \beta, \theta \in \mathbb{R}^+,\label{eqn19}\end{equation}$$

with the Laplace exponent

(4.4) $$\begin{equation}\Phi(u)=-\theta \Gamma (-\alpha) \big[(\beta + u )^\alpha - \beta^\alpha\big],\end{equation}$$

the mean at unit time

$$\begin{equation}\mu_Z=\mathbb{E}[Z_1]=\theta \beta^{\alpha-1} \Gamma(1-\alpha),\end{equation}$$

where $\alpha$ is the stability index, $\theta$ is the intensity parameter, and $\beta$ is the tilting parameter. The tempered stable (TS) subordinator is a Lévy process $\{Z_t\colon t\ge 0 \}$ such that $Z_1 \sim \operatorname{\textsf{TS}} (\alpha,\beta,\theta)$ for $0<\alpha<1$ and $\beta,\theta>0$ .

In fact, tempered stable distribution with three parameters is a very general and flexible distribution. The stable index $\alpha$ determines the importance of small jumps for the process trajectories, the intensity parameter $\theta$ controls the intensity of jumps, and the tilting parameter $\beta$ determines the decay rate of large jumps. In particular, if $\alpha=\frac{1}{2}$ , it reduces to a very important distribution, the inverse Gaussian (IG) distribution, which can be interpreted as the distribution of the first passage time of a Brownian motion to an absorbing barrier. So, this family of tempered stable subordinators in Definition 4.2 also covers the inverse Gaussian (IG) subordinator as an important special case [Reference Barndorff-Nielsen11, Reference Barndorff-Nielsen12]. Conventionally, the IG distribution is denoted by $\operatorname{\textsf{IG}} (\mu_{\operatorname{\textsf{IG}} }, \lambda_{\operatorname{\textsf{IG}} })$ , where $\mu_{\operatorname{\textsf{IG}} }$ is the mean parameter and $\lambda_{\operatorname{\textsf{IG}} }$ is the rate parameter; see [Reference Chhikara and Folks33] for a detailed introduction to IG distributions. The IG subordinator is a special TS subordinator such that the BDLP $Z_t\sim \operatorname{\textsf{IG}} ( \frac{t}{c}, t^2 )$ for any $c,t \in \mathbb{R}^+$ , that is,

$$ \operatorname{\textsf{IG}} \bigg( \dfrac{t}{c},~ t^2\bigg) {\,{\buildrel \mathcal{D} \over =}} \operatorname{\textsf{TS}} \bigg( \dfrac{1}{2}, ~ \dfrac{c^2}{2}, ~ \dfrac{t}{\sqrt{2 \pi} } \bigg). $$

Many scholars have adopted gamma, inverse Gaussian, and tempered stable subordinators as the building blocks to further construct other useful stochastic processes, and there are tremendous relevant papers and work in the literature; see [Reference Barndorff-Nielsen, Shephard, Barndorff-Nielsen, Resnick and Mikosch14, Reference Barndorff-Nielsen and Shephard15, Reference Barndorff-Nielsen and Shephard17], [Reference Cont and Tankov35], [Reference Kyprianou76], [Reference Schoutens and Cariboni95], and [Reference Li and Linetsky81] to name a few. Their main attraction is that the resulting models could possess skewness and leptokurtosis marginally, while remaining highly mathematically tractable. Our new models introduced in this paper can additionally incorporate the ‘contagion’ property, which is also very desirable from the perspective of applications.

Exact simulation of interarrival time. For each case, the interarrival time can be simulated via the general algorithm of thinning scheme, Algorithm 3.2, by simply calculating $\zeta_\infty$ in (3.7) from (3.5) explicitly as

$$\begin{equation}\zeta_\infty=\begin{cases} \displaystyle \varrho a\ln\bigg( 1 +\dfrac{1}{\delta b}\bigg) & \mbox{for {gamma}},\\\\[-5pt] \displaystyle -\varrho\theta\Gamma(-\alpha)\bigg[ \bigg(\beta+\dfrac{1}{\delta}\bigg)^\alpha -\beta^\alpha\bigg] & \mbox{for TS}. \end{cases}\end{equation}$$

Exact simulation of pre-jump intensity. The pre-jump intensity level conditional on the realisation of interarrival time is characterised by the Laplace transforms (3.8) in Theorem 3.2, with Lévy measures $\nu$ and Laplace exponents $\Phi$ specified by (4.1)–(4.2) and (4.3)–(4.4) for the gamma and TS cases, respectively. For both cases, the integral transforms of the pre-jump intensity levels can in fact be broken down into several simple elements.

The Laplace transform of pre-jump intensity level with three terms in (3.8) in fact consists of two parts: (3.11) and (3.12). Strikingly, based on Theorem 3.2, the first two terms of (3.8), i.e. (3.12), can be further exactly decomposed for the specified TS and gamma cases, respectively, as follows.

Algorithm 4.1.

(Exact simulation of pre-jump intensity level for $\Gamma$ -contagion) For the $\Gamma$ -contagion, conditional on the intensity level $\lambda_{T_i}$ and the realisation of the ${(i+1)}$ th interarrival time $\tau_{i+1}=\tau$ , the distribution of the ${(i+1)}$ th pre-jump intensity level $\lambda_{T_i+\tau^-}$ can be exactly decomposed by

(4.5) $$\begin{equation} \lambda_{{T_i}+\tau^-} \!{\mid}\! \lambda_{T_i} {\,{\buildrel \mathcal{D} \over =}\,} w \lambda_{T_i} + \widetilde{\Gamma} + \skew2\widetilde{B}\times S + \sum^{\widetilde{N}}_{j=1} S_j, \label{eqn21} \end{equation}$$

where $\widetilde{\Gamma}, \skew2\widetilde{B}, S, \widetilde{N}$ and $\{S_j\}_{j=1,2,\ldots}$ are all independent of each other,

• $\widetilde{\Gamma}$ is a gamma random variable of

  $$\begin{equation}\widetilde{\Gamma} \sim \operatorname{\textsf{Gamma}}\bigg( -\dfrac{a\varrho }{\delta}\ln w, ~ \dfrac{\vartheta}{w}-\dfrac{1}{\delta} \bigg), \quad \vartheta\,{:\!=}\, b +\dfrac{1}{\delta}{,}\end{equation}$$

• $\skew2\widetilde{B}$ is a Bernoulli random variable taking 0 with probability $p_1$ and 1 with probability $p_2$ , and

  $$\begin{equation}p_1= \dfrac{w\lambda_{T_i}}{\frac{a\varrho }{\delta}C + w\lambda_{T_i}}, \quad p_2 =\dfrac{\frac{a\varrho}{\delta}C}{\frac{a\varrho }{\delta}C + w\lambda_{T_i}}, \quad C \,{:\!=}\, \delta\ln\bigg(\dfrac{b\delta+1-w}{b\delta}\bigg){,}\end{equation}$$

• S is an exponential random variable of $S \sim \operatorname{\textsf{Exp}} ( \vartheta W_0 -\frac{1}{\delta})$ and

  (4.6) $$\begin{equation} W_0 {\,{\buildrel \mathcal{D} \over =}\,} [ 1- b \delta ( {{\rm e}}^{\frac{C}{\delta}U_0} -1 ) ]^{-1}, \quad U_0 \sim \operatorname{\textsf{U}} [0,1]{,} \label{eqn22} \end{equation}$$

• $\widetilde{N}$ is a Poisson random variable of rate $\frac{a\varrho}{\delta}\vartheta C_w$ and

  $$ C_w\,{:\!=}\,\int_1^{\frac{1}{w}}\dfrac{\ln u }{\vartheta u-\frac{1}{\delta}}\,{\rm d}u{,} $$

• $\{S_j\}_{j=1,2,\ldots}$ are i.i.d. with $ S_j \sim \operatorname{\textsf{Exp}} ( \vartheta W -\frac{1}{\delta} )$ , and W can be exactly simulated via the A/R scheme of Algorithm B.1.

Proof.

In fact, Algorithm 4.1 is only an explicit specification of Theorem 3.2. Let us first calculate the first two terms of (3.8), i.e. (3.12), by

$$\begin{equation} \dfrac{\mathbb{E}\big[{{\rm e}}^{-v\lambda_{T_i+\tau^-}}\,{{\rm e}}^{-(\Lambda_{T_i+\tau^-}-\Lambda_{T_i})} \!{\mid}\! \lambda_{T_i}\big]}{\mathbb{E}\big[{{\rm e}}^{-(\Lambda_{T_i+\tau^-}-\Lambda_{T_i})} \!{\mid}\! \lambda_{T_i}\big]} \notag \\[4pt] \qquad = {{\rm e}}^{-vw \lambda_{T_i}} \times \exp\bigg( -\dfrac{a\varrho}{\delta} \ln\bigg(\dfrac{1}{w} \bigg) \int_{0}^{\infty} (1-{{\rm e}}^{-vs} )s^{-1} \,{{\rm e}}^{- (\frac{\vartheta}{w}-\frac{1}{\delta})s}\,{\rm d}s\bigg) \notag \\[4pt] \qquad\quad\, \times\exp\bigg(-\dfrac{a{\vartheta}\varrho}{\delta}\int_0^\infty (1- {{\rm e}}^{-vs} ) \int_1^{\frac{1}{w}}\bigg({\vartheta}u-\dfrac{1}{\delta} \bigg)\,{{\rm e}}^{- ({\vartheta}u-\frac{1}{\delta} )s}\dfrac{\ln u}{ {\vartheta}u-\frac{1}{\delta} }\,{\rm d}u \,{\rm d}s\bigg).\end{equation}$$

Then, by calculating the whole equation (3.8) more explicitly, the conditional Laplace transform of the pre-jump intensity level $\lambda_{T_i+\tau^-}$ can be decomposed into four parts:

(4.7) $$\begin{align}\mathbb{E}\big[ {{\rm e}}^ {-v \lambda_{T_i+\tau_{i+1}^-}} \!{\mid}\! \tau_{i+1}=\tau, \lambda_{T_i}\big] \notag \\[4pt] \qquad = {{\rm e}}^{-vw \lambda_{T_i}} \times \exp\bigg( -\dfrac{\varrho}{\delta}\ln\bigg(\dfrac{1}{w} \bigg)\int_{0}^{\infty} (1-{{\rm e}}^{-vs}) s^{-1} \,{{\rm e}}^{- ( \frac{\theta}{w}-\frac{1}{\delta})s }\,{\rm d}s\bigg) \notag \\[4pt] \qquad\quad\, \times\exp\bigg(-\dfrac{a{\vartheta}\varrho}{\delta}\int_0^\infty (1- {{\rm e}}^{-vs} ) \int_1^{\frac{1}{w}}\bigg({\vartheta}u-\dfrac{1}{\delta} \bigg)\,{{\rm e}}^{- ({\vartheta}u-\frac{1}{\delta} )s}\dfrac{\ln u}{ {\vartheta}u-\frac{1}{\delta} }\,{\rm d}u\,{\rm d}s\bigg) \notag \\[4pt] \qquad\quad\, \times \dfrac{\frac{a\varrho C}{\delta}}{\frac{a\varrho C}{\delta}+w\lambda_{T_i}} \notag \\[4pt] \qquad\quad\,\times \bigg[{\int_{0}^{\infty}\,{{\rm e}}^{-vs}\int_{1}^{\frac{1}{w}}\bigg({\vartheta}u+\dfrac{1}{\delta} \bigg) \,{{\rm e}}^{- ({\vartheta}u-\frac{1}{\delta} )s}\dfrac{1}{C({\vartheta}u^2-\frac{1}{\delta}u )}\,{\rm d}u\,{\rm d}s}+\dfrac{w\lambda_{T_i}}{\frac{a\varrho C}{\delta}+w\lambda_{T_i}}\bigg]. \label{eqn23}\end{align}$$

This decomposition of (4.7) indicates that the conditional distribution of $\lambda_{{T_i}+\tau^-}$ is the sum of four independent simple elements of (4.5): (i) one deterministic trend, (ii) one random variable $\skew2\widetilde{B}\times S$ , (iii) one gamma random variable, and (iv) one compound Poisson random variable. Note that $\skew2\widetilde{B}\times S$ can be alternatively defined as just one single random variable by

$$\begin{equation}\skew2\widetilde{B}\times S{\,{\buildrel \mathcal{D} \over =}\,}\begin{cases}0 & \mbox{with probability } p_1= \dfrac{w\lambda_{T_i}}{\frac{a\varrho }{\delta}C + w\lambda_{T_i}},\\\\[-5pt]S \sim \operatorname{\textsf{Exp}} \bigg( \vartheta W_0 -\dfrac{1}{\delta} \bigg) & \mbox{with probability } p_2 =\dfrac{\frac{a\varrho}{\delta}C}{\frac{a\varrho }{\delta}C + w\lambda_{T_i}}.\end{cases}\end{equation}$$

The CDF of $W_0$ is

$$ F_{W_0}(u) = \dfrac{1}{C_w (a + \frac{1}{\delta})}\ln\bigg( \dfrac{ (a + \frac{1}{\delta} ) u - \frac{1}{\delta} }{a} \bigg), \quad u\in \bigg[1, \dfrac{1}{w} \bigg], $$

which can be inverted explicitly, so we have (4.6). The compound Poisson random variable $\sum_{j=1}^{\tilde{N}} S_j$ has the Laplace transform

$$ \exp\bigg(-\dfrac{a{\vartheta}\varrho}{\delta}C_w \int_0^\infty (1- {{\rm e}}^{-vs} )\int_1^{\frac{1}{w}}\bigg({\vartheta}u-\dfrac{1}{\delta} \bigg)\,{{\rm e}}^{- ({\vartheta}u-\frac{1}{\delta} )s}\dfrac{\ln u }{C_w ({\vartheta}u-\frac{1}{\delta} )} \,{\rm d}u\,{\rm d}s\bigg), $$

so the Poisson rate is $\frac{a\vartheta\varrho}{\delta}C_w$ , and jump-sizes $\{S_j \}_{j=1,2,\ldots}$ follow an exponential distribution with rate $(\vartheta W-\frac{1}{\delta} )$ . Here, W is a well-defined random variable with density (B.1), which can be exactly simulated via the A/R scheme of Algorithm B.1. For more details on designing simulation algorithms based on the A/R mechanism, see [Reference Glasserman65] and Reference Asmussen and Glynn8].

Algorithm 4.2.

(Exact simulation of pre-jump intensity level for TS-contagion.) For the TS-contagion, conditional on the intensity level $\lambda_{T_i}$ and the realisation of the ${(i+1)}$ th interarrival time $\tau_{i+1}=\tau$ , the distribution of the ${(i+1)}$ th pre-jump intensity level $\lambda_{T_i+\tau^-}$ can be exactly decomposed by

(4.8) $$\begin{equation} \lambda_{T_i+ \tau^-} \!{\mid}\! \lambda_{T_i} {\,{\buildrel \mathcal{D} \over =}\,} w \lambda_{T_i} + \widetilde{TS} + \skew2\widetilde{B}\times S + \sum^{\widetilde{N}}_{k=1} S_k, \label{eqn24} \end{equation}$$

where $\widetilde{TS}, \skew2\widetilde{B}, S, \widetilde{N}$ and $\{S_k\}_{k=1,2,\ldots}$ are all independent of each other,

• $\widetilde{TS}$ is an TS random variable of

  (4.9) $$\begin{equation} \widetilde{TS} \sim \operatorname{\textsf{TS}} \bigg( \alpha, ~ \dfrac{\kappa}{w} - \dfrac{1}{\delta}, ~ \dfrac{\theta\varrho } {\alpha\delta} (1- w ^ {\alpha} ) \bigg), \quad \kappa \,{:\!=}\, \beta + \dfrac{1}{\delta}{,} \label{eqn25} \end{equation}$$

• $\skew2\widetilde{B}$ is a Bernoulli random variable taking 0 with probability $p_1$ and 1 with probability $p_2$ , and

  $$\begin{gather*} p_1\,{:\!=}\, \dfrac{w\lambda_{T_i}}{ \frac{\theta \varrho D}{\delta}\Gamma(1-\alpha) + w \lambda_{T_i} }, \quad p_2\,{:\!=}\,\dfrac{\frac{\theta\varrho D}{\delta}\Gamma(1-\alpha)}{ \frac{\theta\varrho D}{\delta} \Gamma(1-\alpha)+ w \lambda_{T_i} }, \\[4pt] D \,{:\!=}\, \dfrac{\delta}{\alpha} \bigg[ \bigg( \kappa-\dfrac{w}{\delta} \bigg)^\alpha - \beta ^ \alpha \bigg]{,} \end{gather*}$$

• S is a mixture-gamma random variable of

  $$ S \sim \operatorname{\textsf{Gamma}}\bigg(1-\alpha, ~ \kappa V_0 - \dfrac{1}{\delta} \bigg) $$
  and
  (4.10) $$\begin{equation} V_0 {\,{\buildrel \mathcal{D} \over =}\,} \bigg[ \delta\beta + 1 - \delta \bigg( \dfrac{\alpha D}{\delta} U_2 + \beta ^{\alpha} \bigg)^{\frac{1}{\alpha}} \bigg]^{-1}, \quad U_2 \sim \operatorname{\textsf{U}} [0,1]{,} \label{eqn26} \end{equation}$$

• $\widetilde{N}$ is a Poisson random variable of rate

  $$ \dfrac{\theta \varrho}{\alpha\delta}\kappa\Gamma(1-\alpha)D_w $$
  and
  $$\begin{equation*} D_w \,{:\!=}\, \int _ 1 ^{\frac{1}{w}} \dfrac{ 1- u ^ {-\alpha} }{( \kappa u -\frac{1}{\delta}) ^ {1-\alpha} } \,{\rm d}u{,} \end{equation*}$$

• $\{S_k\}_{k=1,2,\ldots}$ are i.i.d. with $S_k \sim \operatorname{\textsf{Gamma}}(1-\alpha , ~ \kappa V - \frac{1}{\delta} )$ , and V can be exactly simulated via the A/R scheme of Algorithm C.1.

Proof.

Algorithm (4.2) is another explicit specification of Theorem 3.2. Similarly to the previous gamma case in Algorithm 4.1, given the Lévy measure (4.3), we can identify that (3.11) is the Laplace transform of $\skew2\widetilde{B}\times S$ from the calculation

$$\begin{align*} &\dfrac{ \displaystyle\dfrac{\varrho}{\delta} \int_{0}^{\infty}\,{{\rm e}}^{-vs}\int_{s}^{\frac{s}{w}}\,{{\rm e}}^{\frac{s}{\delta}}\,{{\rm e}}^{-\frac{y}{\delta}}\nu({\rm d}y)\,{\rm d}s+w\lambda_{T_i}} {\displaystyle\dfrac{\varrho}{\delta} \int_{0}^{\infty}\int_{s}^{\frac{s}{w}}\,{{\rm e}}^{\frac{s}{\delta}}\,{{\rm e}}^{-\frac{y}{\delta}}\nu({\rm d}y)\,{\rm d}s+w\lambda_{T_i}}\notag \\ &\qquad = \dfrac{ \displaystyle\dfrac{\theta\varrho}{\delta} \int_0^\infty \,{{\rm e}}^{-vs}\int_1^{\frac{1}{w}}s^{(1-\alpha)-1}\,{{\rm e}}^{- (\kappa u-\frac{1}{\delta} )s}u^{-1-\alpha} \,{\rm d}u\,{\rm d}s + w \lambda_{T_i}} {\displaystyle\int_0^\infty \int_1^{\frac{1}{w}}s^{(1-\alpha)-1}\,{{\rm e}}^{- (\kappa u-\frac{1}{\delta} )s}u^{-1-\alpha} \,{\rm d}u\,{\rm d}s+w\lambda_{T_i}} \notag \\[4pt] &\qquad = \dfrac{w\lambda_{T_i}}{ \frac{\theta\varrho D}{\delta}\Gamma(1-\alpha)+w\lambda_{T_i}} \times 1 + \dfrac{ \frac{\theta\varrho D}{\delta}\Gamma(1-\alpha)} { \frac{\theta\varrho D}{\delta}\Gamma(1-\alpha)+w\lambda_{T_i}}\\[4pt] & \qquad\quad\, \times \int_0^\infty \,{{\rm e}}^{-vs}\int_1^{\frac{1}{w}}\dfrac{(\kappa u-\frac{1}{\delta})^{1-\alpha}}{\Gamma(1-\alpha)}s^{(1-\alpha)-1} \,{{\rm e}}^{- (\kappa u-\frac{1}{\delta})s}\dfrac{u^{-1-\alpha}}{D (\kappa u-\frac{1}{\delta})^{1-\alpha}}\,{\rm d}u\,{\rm d}s \notag \\&\qquad = p_1 \times \mathbb{E}[ {{\rm e}}^ {-v 0} ] + p_2 \times \mathbb{E}[ {{\rm e}}^ {- v S} ], \notag \end{align*}$$

where

$$D=\int _ 1 ^ {\frac{1}{w} } \dfrac{u^{-1-\alpha}}{(\kappa u-\frac{1}{\delta})^{1-\alpha}}\,{\rm d}u=\dfrac{\delta}{\alpha} \bigg[ \bigg( \kappa-\dfrac{w}{\delta} \bigg)^\alpha - \beta ^ \alpha \bigg].$$

Hence, the outcome of $\skew2\widetilde{B}\times S$ is trivially equal to 0 with probability $p_1$ or the random variable S with probability $p_2$ . S follows a mixture-gamma distribution with shape parameter $1-\alpha$ and rate parameter $\kappa V_0 - \frac{1}{\delta}$ . Here, $V_0$ is a well-defined random variable with density function

$$\begin{equation}f_{V_0}(u)=\dfrac{ u ^ {-\alpha -1} }{D ( \kappa u -\frac{1}{\delta} )^ {1-\alpha} },\quad u\in \bigg[1, \dfrac{1}{w} \bigg].\end{equation}$$

It can be directly simulated via the explicit inverse transform (4.10), as its CDF is

$$ F_{V_0}(u) = \dfrac{\delta}{\alpha D} \bigg[ u ^ {-\alpha} \bigg( \kappa u -\dfrac{1}{\delta}\bigg)^\alpha -\beta ^ \alpha\bigg], \quad u\in \bigg[1, \dfrac{1}{w} \bigg]. $$

The first two terms of (3.8), i.e. (3.12), can be expressed by

$$\begin{equation}\dfrac{\mathbb{E}\big[{{\rm e}}^{-v\lambda_{T_i+ \tau^-}}\,{{\rm e}}^{-(\Lambda_{T_i+ \tau^-}-\Lambda_{T_i})} \!{\mid}\! \lambda_{T_i}\big]}{\mathbb{E}\big[{{\rm e}}^{-(\Lambda_{T_i+ \tau^-}-\Lambda_{T_i})} \!{\mid}\! \lambda_{T_i}\big]} \notag \\[4pt] \qquad = {{\rm e}}^{-vw\lambda_{T_i}} \times \exp\bigg(-\dfrac{\theta\varrho}{\alpha\delta} (1-w^\alpha ) \int_{0}^{\infty}(1-{{\rm e}}^{-vs} )\dfrac{1}{s^{\alpha+1}} \,{{\rm e}}^{- (\frac{\kappa}{w}-\frac{1}{\delta})s}\,{\rm d}s\bigg) \notag \\[4pt] \qquad\quad\, \times \exp\bigg(-\dfrac{\theta\varrho}{\alpha\delta}\kappa\Gamma(1-\alpha) \,{I} \bigg),\end{equation}$$

where

$$ I = \int_{0}^{\infty} (1-{{\rm e}}^{-vs} )\int_1^{\frac{1}{w}}\dfrac{(\kappa u-\frac{1}{\delta})^{1-\alpha}}{\Gamma(1-\alpha)}s^{(1-\alpha)-1}\dfrac{1-u^{-\alpha}}{(\kappa u-\frac{1}{\delta})^{1-\alpha}}\,{\rm d}u\,{\rm d}s , $$

since the Lévy measure $\nu$ for the TS subordinator is specified in (4.3), and the Laplace exponent of (3.12) can be rewritten as

$$\begin{align*} & \varrho\int_{v}^{\frac{1}{\delta}+(v-\frac{1}{\delta})w}\dfrac{\Phi \big(u-\tfrac{u-\frac{1}{\delta}}{v-\frac{1}{\delta}}v\big)-\Phi(u)}{1-\delta u}\,{\rm d}u \notag \\[4pt] &\qquad = \dfrac{\varrho}{\delta}\int_{0}^{\infty} (1-{{\rm e}}^{-vs})\dfrac{1}{s}\,{{\rm e}}^{\frac{s}{\delta}} \int_{s}^{\frac{s}{w}}\dfrac{\theta}{y^{\alpha+1}}\,{{\rm e}}^{-\kappa \frac{s}{w}}\,{\rm d}y\,{\rm d}s \notag \\[4pt] &\qquad\quad\, + \dfrac{\varrho}{\delta}\int_{0}^{\infty} (1-{{\rm e}}^{-vs})\dfrac{ {{\rm e}}^{\frac{s}{\delta}}}{s} \int_{s}^{\frac{s}{w}}\dfrac{\theta ( {{\rm e}}^{-\kappa y}-{{\rm e}}^{- \frac{\kappa s}{w}})}{y^{\alpha+1}} \,{\rm d}y\,{\rm d}s \notag \\[4pt] &\qquad = \int_{0}^{\infty} (1-{{\rm e}}^{-vs})\dfrac{\theta\varrho}{\alpha\delta} (1-w^\alpha ) \dfrac{1}{s^{\alpha+1}}\,{{\rm e}}^{- (\frac{\kappa}{w}-\frac{1}{\delta})s}\,{\rm d}s \notag \\[4pt] &\qquad\quad\, +\dfrac{\theta\varrho}{\alpha \delta}\kappa\Gamma(1-\alpha) \,{J} \end{align*}$$

where

$$ J=\int_0^\infty (1 -{{\rm e}}^{-vs} ) \int_1^{\frac{1}{w}}\dfrac{ (\kappa u-\frac{1}{\delta})^{1-\alpha}}{\Gamma(1-\alpha)}s^{(1-\alpha)-1}\,{{\rm e}}^{- (\kappa u-\frac{1}{\delta})s}\dfrac{1-u^{-\alpha}}{(\kappa u-\frac{1}{\delta})^{1-\alpha}}\,{\rm d}u\,{\rm d}s . $$

So, for (3.12), it consists of three components: one deterministic trend, a TS process, and a compound Poisson process. In particular, the rate of the compound Poisson process is $\widetilde{N}$ is

$$ \dfrac{\theta\varrho}{\alpha\delta}\kappa\Gamma(1-\alpha)D_w, $$

and the jump sizes follow a mixture-gamma distribution with shape parameter $(1-\alpha)$ and rate parameter $\kappa V -\frac{1}{\delta}$ . Here, V is a well-defined random variable with density

(4.11) $$\begin{equation} f_{V}(u) = \dfrac{1-u ^ {-\alpha } }{D_w (\kappa u-\frac{1}{\delta} ) ^ {1-\alpha} }, \quad u \in \bigg[1,\dfrac{1}{w} \bigg]. \label{eqn27} \end{equation}$$

Overall, we have the conditional Laplace transform of pre-jump intensity level explicitly as

$$\begin{equation}\mathbb{E}\big[ {{\rm e}}^ {-v \lambda_{T_i+\tau_{i+1}^-}} \!{\mid}\! \tau_{i+1}=\tau, \lambda_{T_i}\big] \notag \\[4pt] \qquad = {{\rm e}}^{-vw\lambda_{T_i}} \times \exp\bigg(-\dfrac{\theta\varrho}{\alpha\delta} (1-w^\alpha ) \int_{0}^{\infty} (1-{{\rm e}}^{-vs} )\dfrac{1}{s^{\alpha+1}} \,{{\rm e}}^{- (\frac{\kappa}{w}-\frac{1}{\delta})s}\,{\rm d}s\bigg) \notag \\[4pt] \qquad\quad\, \times \exp\bigg(-\dfrac{\theta\varrho}{\alpha\delta}\kappa\Gamma(1-\alpha)D_w \,{K}\bigg) \\[4pt] \qquad\quad\, \times \bigg[\dfrac{w\lambda_{T_i}}{\frac{\theta\varrho D}{\delta}\Gamma(1-\alpha)+w\lambda_{T_i}}+\dfrac{\frac{\theta\varrho D}{\delta}\Gamma(1-\alpha)} {\frac{\theta\varrho D}{\delta}\Gamma(1-\alpha)+w\lambda_{T_i}} \notag \\[4pt] \qquad\quad\, \times \int_0^\infty \,{{\rm e}}^{-vs}\int_1^{\frac{1}{w}}\dfrac{(\kappa u-\frac{1}{\delta})^{1-\alpha}}{\Gamma(1-\alpha)}s^{(1-\alpha)-1} \,{{\rm e}}^{- (\kappa u-\frac{1}{\delta} )s}\dfrac{u^{-1-\alpha}}{D (\kappa u-\frac{1}{\delta})^{1-\alpha}}\,{\rm d}u\,{\rm d}s \bigg],\end{equation}$$

where

$$ K =\int_{0}^{\infty} (1-{{\rm e}}^{-vs} )\int_1^{\frac{1}{w}}\dfrac{(\kappa u-\frac{1}{\delta})^{1-\alpha}}{\Gamma(1-\alpha)}s^{(1-\alpha)-1}\dfrac{1-u^{-\alpha}}{D_w (\kappa u-\frac{1}{\delta})^{1-\alpha}}\,{\rm d}u\,{\rm d}s. $$

We can identify from the Laplace transforms above that the distribution of the ${(i+1)}$ th pre-jump intensity level $\lambda_{T_i + \tau^-}$ conditional on $\lambda_{T_i}$ is exactly equal in distribution to the sum of four simple elements provided in (4.8). All these components can be simulated exactly. To simulate the TS random variable, one could use existing algorithms provided in [Reference Brix23], [Reference Devroye53], [Reference Hofert73] or backward recursive (BR) scheme provided in [Reference Dassios, Qu and Zhao50]. And to sample the compound Poisson random variable $\widetilde{N}$ , one first needs to generate the intermediate random variable V with density (4.11). Since there is no closed form for the inverse function of the CDF of V, we have to rely on the A/R scheme of Algorithm C.1.

Exact simulation of self-exciting jumps. Conditional on the realisations of the interarrival time and pre-jump intensity level as above, the associated self-exciting jump can easily be simulated by just following Algorithm 3.3 in general both for gamma and TS contagion processes.

5. Numerical experiments

In this section, let us take the TS-contagion model as an example to illustrate the performance of our exact scheme via extensive numerical experiments, and postpone the implementation of the gamma-contagion model to Section 6 with more financial applications. The simulation experiments here and elsewhere in this paper are all conducted on a desktop PC with an Intel^® Core^TM i7-3770S CPU@3.10 GHz processor, 8.00 GB RAM, Windows 7^®, 64-bit Operating System; the algorithms are coded and performed in MATLAB^® (R2014a), and the computation time is measured by the elapsed CPU time in seconds. The true value of the conditional expectation of $N_T$ for any fixed time $T>0$ provided in Proposition 2.1 is used to numerically validate and test our algorithms. The associated errors from the true values are reported by three standard measures.

1. (i) $Error= \text{estimated value} - \text{true value}$ .

2. (ii) $\textit{Relative error}\ (\textit{error}\ \%) = \dfrac{\text{estimated value} - \text{true value}}{\text{true value}}$ .

3. (iii) Root mean square error $\text{RMSE}=\sqrt{\text{bias}^2 + \text{SE}^2}$ , where the SE is the standard error of simulation output, and the bias is the difference between the expectation of the estimator and the associated true (theoretical) value. For the algorithm of exact simulation here, the bias is conventionally set to zero.

We implement Algorithm 3.3 for the TS and IG cases in a fixed period of [0, T] with and without self-exciting jumps.

• Case I. Jump process with non-Gaussian OU intensity (see Definition 2.1).

• Case II. Self-exciting jump process with non-Gaussian OU intensity (of Definition 2.2).

Note that as an intermediate step we have to generate the random variable $\widetilde{TS}$ of (4.9) in Algorithm 4.2 for the general TS case. Several algorithms are available in the literature, including approximation-based algorithms such as infinite series representation [Reference Rosiski, Barndorff-Nielsen, Resnick and Mikosch91] and exact algorithms such as simple stable rejection (SSR) [Reference Brix23], the double rejection algorithm [Reference Devroye53], the fast rejection algorithm [Reference Hofert73], and the backward recursive (BR) scheme [Reference Dassios, Qu and Zhao50]. However, the choice of algorithms is in fact not the focus of our paper. For the purpose of demonstration, we adopt the BR scheme. It works extremely efficiently for some families, including those with stability index the binary fractions $\frac{1}{2^n}$ , $n=1,2,\ldots,$ which can easily be simulated by recursively generating IG random variables without the A/R mechanism. For a general parameter setting, one could simply adopt other algorithms such as the SSR scheme of Algorithm E.1.

For numerical implementation, we further assume that the sizes of self-exciting jumps follow an exponential distribution of rate $\gamma>0$ , i.e. $X_i \sim \operatorname{\textsf{Exp}} (\gamma)$ , and the stable index takes the values of, say, $\alpha=1/4$ for the TS case and $\alpha=1/2$ for the IG case, respectively. $\widetilde{TS}$ of (4.9) is simulated using the BR scheme specifically designed for $\alpha=1/4$ in Algorithm D.1, and the parameters are set as follows.

• TS Case I. $(\delta, \varrho; \alpha, \beta, \theta; \lambda_0)=(1.0,0.5;0.25,0.2,0.25; 0.5)$ .

• TS Case II. $(\delta, \varrho; \alpha, \beta, \theta; \gamma; \lambda_0)=(1.0,0.5;0.25,0.2,0.25; 5.0; 0.5)$ .

• IG Case I. $(\delta, \varrho; c; \lambda_0)=(1.0,0.5;0.5; 0.5)$ .

• IG Case II. $(\delta, \varrho; c; \gamma; \lambda_0)=(1.0,0.5;0.5; 4.0; 0.5)$ .

Simulated sample paths of the point processes within a long period of $t\in[0, 500]$ with the associated histograms are plotted in Figure 3, where the clustering or ‘contagious’ arrivals of jumps can be clearly presented. Furthermore, to measure the accuracy and efficiency of our scheme, we carry out the convergence analysis: Figure 4 presents log-log plots for the RMSE against the CPU time for each case in two different time horizons $T=2, 5$ , respectively, and the associated results in detail are reported in Tables 1 and 2. Overall, from the numerical results reported in this section, it is evident that our exact scheme can achieve a very high level of accuracy and efficiency.

Fig. 3.

Simulated sample paths of the point processes and the associated time-series plots for Cases I and II: TS Case I $(\delta, \varrho; \alpha, \beta, \theta; \lambda_0)=(1.0,0.5;0.25,0.2,0.25; 0.5)$ , TS Case II $(\delta, \varrho; \alpha, \beta, \theta; \gamma; \lambda_0)=(1.0,0.5;0.25,0.2,0.25; 5.0; 0.5)$ ; IG Case I $(\delta, \varrho; c; \lambda_0)=(1.0,0.5;0.5; 0.5)$ , IG Case II $(\delta, \varrho; c; \gamma; \lambda_0)=(1.0,0.5;0.5; 4.0; 0.5)$ .

Fig. 4.

Convergence analysis via RMSE versus CPU time by log-log plots for Cases I and II: TS Case I $(\delta, \varrho; \alpha, \beta, \theta; \lambda_0)=(1.0,0.5;0.25,0.2,0.25; 0.5)$ , TS Case II $(\delta, \varrho; \alpha, \beta, \theta; \gamma; \lambda_0)=(1.0,0.5; 0.25,0.2,0.25;5.0; 0.5)$ ; IG Case I $(\delta, \varrho; c; \lambda_0)=(1.0,0.5;0.5; 0.5)$ , IG Case II $(\delta, \varrho; c; \gamma; \lambda_0)=(1.0,0.5; 0.5; 4.0; 0.5)$ .

Table 1.

Simulation results for Cases I and II: TS Case I $(\delta, \varrho; \alpha, \beta, \theta; \lambda_0)=(1.0,0.5; 0.25,0.2,0.25; 0.5)$ TS Case II $(\delta, \varrho; \alpha, \beta, \theta; \gamma; \lambda_0)=(1.0,0.5;0.25,0.2,0.25; 5.0; 0.5)$ .

Table 2.

Simulation results for Cases I and II: IG Case I $(\delta, \varrho; c; \lambda_0)=(1.0,0.5;0.5; 0.5)$ , IG Case II $(\delta, \varrho; c; \gamma; \lambda_0)=(1.0,0.5;0.5; 4.0; 0.5)$

To be even more prudent, the simulation for the interarrival time based on the simplified thinning scheme of Algorithm 3.2, as an intermediate step, can also be tested separately. To numerically assess its accuracy and efficiency, we compare the simulated results of V ^* with its theoretical tail distribution $\mathbb{P } \{ V^* >\tau \}$ as specified in (3.6), which can be calculated explicitly by substituting the Laplace exponent $\Phi$ from (4.4). We set the parameters by $(\delta, \varrho; \alpha, \beta, \theta)=(0.5,1.0;0.9,0.2,0.25)$ , and each estimation is based on 10⁵ replications. Error percentages (Error%) for measuring relative errors are reported in Table 3. The total CPU time for producing the whole Table 3 is only 12.64 seconds, and the error percentages are all very tiny.

Table 3.

Comparison between the theoretical formulas and the associated simulation results for the simplified thinning scheme of Algorithm 3.2 with each estimation based on 10⁵ replications.

6. Comprehensive risk analysis for a large portfolio facing contagious losses and unexpected exogenous gamma shocks

It has now been widely recognised among academics and financial practitioners that risk spreads through highly interconnected business networks, and defaults could trigger more defaults via a ‘domino’ effect. The resulting losses presented in financial markets could be amplified. As earlier mentioned in Section 4, gamma distribution is a popular building block in financial applications. In particular, it plays an important role in credit risk modelling. For instance, both the widely used framework of [42] in the banking industry and influential papers by [Reference Gordy67], [Reference Gordy68, and Reference Elsinger, Lehar and Summer59] in the literature, assumed that macroeconomic factors are driven by independent gamma-distributed random variables. Further, it can also be equipped as a fundamental risk driver for price movements, e.g. the popular variance gamma model [Reference Madan and Seneta83, Reference Madan, Carr and Chang84]. More recent evidence has been found by [Reference Giesecke, Longstaff, Schaefer and Strebulaev64] that long periods with relatively few defaults are followed by episodes of significant clustering of defaults, and the resulting distribution of default rates is highly skewed towards large values. This motivates us to adopt our new model of gamma contagion as an example for applications in risk management for a portfolio facing a ‘domino’ effect of losses. We assume that exogenous commonly shared risk is dynamically powered by a gamma process. More precisely, we adopt the gamma distribution as the fundamental driver of randomness (or gamma shock) to construct the OU- $\Gamma$ interarrival intensity for a point process $N_t$ in Definition 2.2. A simulated path of this interarrival intensity process within the time period $t \in[0, 5]$ based on the parameter setting $ (\delta, \varrho; a, b; \lambda_0 )=(1.0,1.0;4.0,0.5; 2.0)$ was plotted in Figure 1 above, where we can observe relatively high-frequency and small shocks. In fact, Barndorff-Nielsen and Shephard [Reference Barndorff-Nielsen and Shephard17] called it the OU- $\Gamma$ process, and it has become a very popular tool for modelling stochastic volatilities in a continuous-time set-up. Hainaut and Devolder [Reference Hainaut and Devolder69] used it as a special case of Cox processes to model human mortality rates, and applied to actuarial valuation in insurance. Eberlein, Madan, Pistorius, and Yor [Reference Eberlein, Madan, Pistorius and Yor57] treated it as a one-factor model for describing the evolution of instantaneous interest rates.n reality, contagion may be triggered by losses or defaults of banks or other financial institutions through inter-institutional lendings in the interbank market, or it may be further amplified due to some common asset holdings of overlapping portfolios [Reference Caccioli, Shrestha, Moore and Farmer27]. We offer some numerical examples of comprehensive risk analysis for a large portfolio facing contagious defaults and losses. We construct a simple contagious loss process to capture the propagated defaults for a generic large pool of financial institutions (banks for short) within a financial system. The following framework, of course, would work generically for other similar institutions. The aggregate loss process of this large portfolio by time t is

$$ L_t = \sum^{N_t}_{i=1} L_i, \quad t\geq 0, $$

where $N_t$ is a $\Gamma$ -contagion process, and $L_i \geq 0$ is the absolute value of the loss size for the ith default, of which the mean is denoted by $\mu_L\,{:\!=}\,\mathbb{E}[L_i]$ for any index i. We assume that the sizes of self-exciting jumps in (2.3) generally satisfy

(6.1) $$\begin{equation}X_i = \varpi_i \times g(H_i),\end{equation}$$

where

• $H_i$ is the history of the loss path until time t, i.e. $H_i\,{:\!=}\, \{L_j\}_{j< i} \cup \{T_j\}_{j\leq i}$ ,

• $\varpi_i>0$ is the amplification multiplier, the amplification mechanisms during financial crises in detail were described and analysed in [Reference Brunnermeier25], [Reference Brunnermeier and Pedersen26], and more recently in [Reference Ahnert and Kakhbod2], which might be dependent on the degree of financial connectivity of the underlying company i to others, or the effects of policymakers’ interventions to limit the extent of contagion,

• $g({\cdot})$ is a general non-negative function of the losses. Since the jump sizes in (6.1) follow a $\mathcal{F}_{T^{-}_i}$ -measurable distribution as explained earlier in Definition 2.2, similarly the function $g({\cdot})$ is also very general, and it can include all past losses, that is,

  $$ g({\cdot}) = g ( T_1, T_2,\ldots,T_i, L_1,\ldots,L_{i-1} ), $$
  which could adopt functional forms, for example, assigning weight to each loss size according to its loss time or similarly to autoregressive time series models.

In fact, (6.1) provides a channel for contagion (or feedback) effects of market participants’ reactions to adverse scenarios. The economic interpretation for this model is that the impacts and the timing of unexpected exogenous gamma shocks acting on the entire portfolio as macroeconomic scenarios are modelled by a mean-reverting gamma-driven OU process. The shocks might not lead to an immediate default but act on the underlying intensity via a positive jump, which increases the default probability afterwards. Meanwhile, endogenic shocks, i.e. contagious losses due to the propagated defaults, are modelled by self-exciting jumps, and the associated magnitudes can be captured by jump sizes $ \{ X_i \}_{i=1,2,\ldots}$ .

The great flexibility of our exact simulation scheme allows us to accurately and efficiently generate highly comprehensive scenarios for risk assessment. In general, our algorithms can simulate sample paths when loss sizes $L_i$ may depend on the entire history of $N_t$ and $\lambda_t$ before or at time $T_i$ . We discuss several circumstances that can be captured by our models as follows.

6.1. A simple benchmark model

The loss occurring within a financial institution may spread via various business channels and eventually trigger subsequent losses of others in markets. Intuitively, a larger loss may make a larger impact. For convenience of illustration, we assume that the sizes of self-exciting jumps satisfy

$$ X_i=\bar{\varpi}\times L_i, $$

where $\bar{\varpi}>0$ is the average amplification multiplier, meaning that each investment has a linear and homogeneous amplification effect. We further assume that each loss size is exponentially distributed, i.e. $L_i\sim \operatorname{\textsf{Exp}} (\ell)$ , $\ell>0$ with mean $\mu_{L}\,{:\!=}\,1/\ell$ . To assess the overall risk of this portfolio, we implement the exact simulation of Algorithms 3 and 4 with (Case I) and without (Case II) contagion in the fixed time period [0, t], respectively:

Case I. $(\delta, \varrho; a, b; \lambda_0)=(0.5,0.5;0.5,2.0; 0.5)$ ,

Case II. $(\delta, \varrho; a, b; \ell, \bar{\varpi}; \lambda_0)=(0.5,0.5;0.5,2.0; 8.0, 2.0; 0.5)$ .

We concentrate on the default number $N_t$ in the system. Case I or II can be considered as a benchmark model, as by Proposition 2.1 the expected default number has analytical forms.

Proposition 6.1.

[Expectation of $N_t$ ] The expected default number until time t is given by

(6.2) $$\begin{equation} \mathbb{E}[N_t \!{\mid}\! \lambda_0] = \lambda_0\dfrac{ 1 - {{\rm e}}^{-\eta t} }{\eta} + \dfrac{\varrho}{\eta}\bigg(t - \dfrac{ 1 - {{\rm e}}^{-\eta t} }{\eta} \bigg)\dfrac{a}{b}, \quad \eta \neq 0, \label{eqn29} \end{equation}$$

where

$$\begin{equation*} \eta= \begin{cases} \delta \mbox{for Case I},\\ \delta-\bar{\varpi}/\ell \mbox{for Case II}. \end{cases} \end{equation*}$$

To explore the models, let us first carry out a sensitivity analysis for the expected default number $\mathbb{E}[N_t \!{\mid}\! \lambda_0]$ with respect to their key parameters for controlling the external gamma shocks, a and b, with and without contagion, and the results are provided respectively in Figure 5. Numerical tests for our algorithms are based on the true means (6.2). The associated errors reported by three standard measures are reported in Table 4. Convergence analysis via log-log plots of the RMSE against the CPU time for Cases I and II and $t=2,5$ is presented in Figure 6. We can observe that simulations are pretty fast with very tiny errors, which provides numerical evidence of accuracy and efficiency for our algorithms.

Fig. 5.

Sensitivity analysis for the expected number of defaults $\mathbb{E}[N_{t=5} \!{\mid}\! \lambda_0]$ with respect to a and b based on the parameters $(\delta, \varrho; \bar{\varpi}; \lambda_0)=(0.5,0.5; 2.0; 0.5)$ .

Table 4.

Simulation results for Cases I and II and time $t=2,5$ : $(\delta, \varrho; a, b; \lambda_0)=(0.5,0.5;0.5,2.0; 0.5)$ for Case I, and $(\delta, \varrho; a, b; \ell, \bar{\varpi}; \lambda_0)=(0.5,0.5;0.5,2.0; 8.0, 2.0; 0.5)$ for Case II

Fig. 6.

Convergence analysis via the log-log plots of the RMSE versus the CPU time for Cases I and II and $t=2,5$ : $(\delta, \varrho; a, b; \lambda_0)=(0.5,0.5;0.5,2.0; 0.5)$ for Case I, and $(\delta, \varrho; a, b; \ell, \bar{\varpi}; \lambda_0)=(0.5,0.5;0.5,2.0; 8.0, 2.0; 0.5)$ for Case II.

6.2. A model with contagion threshold

In reality, each loss might not necessarily cause contagion immediately throughout the entire system. Contagion may be only triggered when the loss exceeds a certain high level, that is, contagion is likely to occur only in severe scenarios, which has also been reported in [Reference Elsinger, Lehar and Summer59]. This circumstance could be modelled by a mixture of Cases I and II, by assuming that the sizes of self-exciting jumps $X_i$ satisfy

$$X_i= \varpi_i \times (L_i-K_i)^+,$$

where $K_i\geq 0$ is the contagion threshold (i.e. the threshold that triggers the contagion effect of the ith loss $L_i$ ), and the contagion has been partially capped. If a bank is more vulnerable, its threshold is more easily reached. Alternatively, we may interpret $K_i$ as a capital buffer, and it could be a certain quantile of the loss distribution $L_i$ . If we assign the same quantile to all banks, it is equivalent to an identical economic capital applying to all banks, which is the assumption made by [Reference Elsinger, Lehar and Summer59]. If the magnitude of loss overshoots the threshold, the bank may become insolvent, and this risk may then spread to other banks (through the interbank market) resulting in an increase in the default intensity of the entire system (but this would not immediately cause other defaults). If the thresholds are very high compared to the levels of loss, then it corresponds to a ‘weak contagion’ environment, whereas if the thresholds are very low, then it is a ‘strong contagion’ environment.

With the contagion threshold, contagion could be partially or fully triggered. Here, for numerical illustration, we assume that losses are exponentially distributed and the amplification multipliers and contagion thresholds are homogeneous, i.e. $L_i \sim \operatorname{\textsf{Exp}} (\ell)$ , $a_i \equiv \bar{\varpi}$ , and $K_i\,{\equiv}\, K \geq0$ . The expected default number can hardly capture the full picture of the risk, and we have to look at the entire distribution. We choose $K=\infty, 1/8, 0$ and plot the estimated probability mass function (PMF) of the total default number within the period of [0, t] in Figure 7, and the corresponding quantiles are reported in Table 5. More specifically, cases $K\,{=}\,\infty, 1/8, 0$ correspond to the non-contagion (i.e. Case I), partial contagion and full contagion (i.e. Case II), respectively. We can clearly observe that when K decreases, the contagion would become more pronounced and the tail of losses become heavier. The system could be more susceptible to contagion risk when capital buffer K is eroded, and contagion effects magnify the content of risk. As summarised by [Reference Eisenberg and Noe58], bank defaults may be driven by losses from market and credit risk (i.e. fundamental default), and bank defaults may, however, also be initiated by contagion as a consequence of other bank failures in the system (i.e. contagious default). The two types of default under our contagion model (i.e. the self-exciting jump sizes are not all equal to zero) are in fact mixed, and interact.

Fig. 7.

Probability mass function (PMF) of the default number $N_{t=5}$ , estimated from 10⁶ replications based on the parameter setting $(\delta, \varrho; a, b; \ell, \bar{\varpi}; \lambda_0)=(0.5,0.5;0.5,2.0; 8.0, 2.0; 0.5)$ , with homogeneous contagion thresholds $K=\infty, 1/8, 0$ , respectively. The associated quantiles are reported in Table 5.

Table 5.

Quantiles of the default number $N_{t=5}$ , estimated from 10⁶ replications based on the parameter setting $(\delta, \varrho; a, b; \ell, \bar{\varpi}; \lambda_0)=(0.5,0.5;0.5,2.0; 8.0, 2.0; 0.5)$ , with homogeneous contagion thresholds $K\,{=}\,\infty, 1/8, 0$ , respectively

6.3. A model with explosive defaults

Contagion or feedback effects could be even further reinforced due to highly leveraged positions (e.g. complicated credit derivatives), and the resulting system thereby becomes explosive; see [Reference Corsi, Marmi and Lillo38] for discussion of the impacts of financial innovations. This scenario would be extremely severe, rare but possible, that is, the entire system is unstable and near the crash boundary. Mathematically, it corresponds to the non-stationary case when $\eta<0$ in our models. This may be due to ‘liquidity black holes’ or ‘fire sales’ of assets: these further depress prices and lead to a sharp drop in liquidity and may also cause other institutions to fail in a self-reinforcing vicious spiral [Reference Cont and Wagalath36,Reference Cont and Wagalath37,Reference Krishnamurthy75,Reference Morris and Shin86]. All previous examples were conducted under the stationary condition $\eta>0$ , and in fact our algorithms can also deal with non-stationary cases. In Figure 8, we offer three representative examples of $\eta=1/4$ , 0, and -1/2 (or $\ell=8$ , 4, and 2) for stationary, critical, and explosive phases, respectively, and the associated quantiles are reported in Table 6. In particular, $\eta=0$ is the critical level of stability. The resulting loss distributions could present heavy tails, which might be very desirable for many regulators and practitioners.

Fig. 8.

Probability mass function (PMF) of the default number $N_{t=5}$ , estimated from 10⁶ replications based on the parameter setting $(\delta, \varrho; a, b; \bar{\varpi}; \lambda_0)=(0.5,0.5;0.5,2.0; 2.0; 0.5)$ , with $\ell=8, 4, 2$ , respectively. The associated quantiles are reported in Table 6.

Table 6.

Quantiles of the default number $N_{t=5}$ , estimated from 10⁶ replications based on the parameter setting $(\delta, \varrho; a, b; \bar{\varpi}; \lambda_0)=(0.5,0.5; 0.5,2.0; 2.0; 0.5)$ , with $\ell=8, 4, 2$ , respectively

7. Conclusion

In this paper we have introduced a new family of self-exciting jump processes whose intensities are driven by non-Gaussian OU processes, namely Lévy-driven contagion processes. Backed by the very large family of Lévy subordinators, it indeed offers much richer choices beyond the classical Hawkes process for modelling the ‘contagion’ of event arrivals in a continuous-time set-up in finance, economics, and many other fields. We have derived some important distributional properties of these new processes which lead to an exact simulation framework in general. In particular, we have developed exact simulation algorithms by the decomposition approach for the gamma and tempered stable cases as typical examples. The algorithms are accurate and efficient, and have been numerically verified and tested by extensive numerical experiments. We also provide applications to portfolio risk management, which again illustrate the efficiency, accuracy, applicability, and flexibility of our algorithms. As a class of reduced-form models, it could easily be extended to pricing financial derivatives, particularly multiple-name credit products (e.g. collateralised debt obligations and mortgage-backed securities). It can be employed empirically when input data are available for parameter calibration. Furthermore, it could be widely applied to many other areas, for example to describe high-frequency trading data in market microstructure, claim arrivals for an insurance portfolio, or jump propagation and disclosure dynamics in financial markets [Reference Acharya, DeMarzo and Kremer1]. Their statistical inference and econometric analysis for this new framework, and further extensions to multi-dimensional point processes for modelling multilateral contagion, as well as further applications and empirical work for portfolio credit risk analysis, could be very interesting and meaningful topics for future research.