The digitalization of business and economic activities has significantly increased the frequency and impact of cyber events around the globe, with alarming consequences to our interconnected world. Emerging cyber events bring about catastrophic and far-reaching losses due to the business interruption. According to the cybersecurity ventures official report (Morgan, 2022), global cybercrime costs could reach US10.5trillion by 2025.

Interconnected business and economic activities have also increased systemic risk events, and one of its examples is the collapse of Lehman Brothers Holdings Inc. in 2008. To avoid economic spillovers between countries and financial sector, the U.S. government intervened in March 2023 to guarantee deposit accounts at Silicon Valley Bank, aiming to restore confidence in the global banking sector. In 2008, the near-failure of the world insurance giant American international group due to the selling of credit default swaps – financial instruments that act like insurance contracts on bonds–could have led to a systemic collapse of the global financial system if not for government assistance. To address interconnections and spillovers between the risks involved and their contagious loss aggregation, a multivariate process is required. For that purpose, we introduce a multivariate compound dynamic contagion process (MCDCP) in this paper.

The dynamic contagion process (DCP) introduced by Dassios and Zhao (2011) is a generalization of the externally exciting Cox process with shot noise intensity and the self-exciting Hawkes process. They analyzed its theoretical distributional properties and applied to credit risk. Compound modeling with univariate self-exciting Hawkes processes can be noticed in Dassios and Zhao (2017), where they developed algorithms for a generalized self-exciting point process with CIR-type intensities, presenting its first moment. Gao et al. (2018) applied the joint Laplace transform of the classical Hawkes process and the compound process in dark pool trading, which does not display bid and ask quotes to the public, presenting its first and second moment.

Bivariate modeling with self-exciting Hawkes processes can be also noticed in Jang and Dassios (2013), where they introduced a bivariate shot noise self-exciting process that can be used for the modeling of catastrophic losses. The joint moment, covariance, and correlation coefficient on the bivariate shot noise self-exciting process can be found in this paper. A multivariate extension on the DCP made by Dong (2014), where the bivariate DCP was introduced, including the cross-exciting contagion effect, and its exact simulation algorithm was provided. Bessy-Roland et al. (2021) proposed a multivariate Hawkes framework for modeling and predicting the frequency of cyber-attacks.

Jang and Oh (2021) studied the compound dynamic contagion process (CDCP) and provided analytic expressions of the Laplace transform/probability generating function of this process. However, to date, no work has been done on its multivariate process. This study focuses on a multivariate compound dynamic contagion process (MCDCP) to quantify the aggregate losses arising from contagious catastrophic events. It provides steps for analyzing the data, with further iterations of the model on real-life datasets being an area for future research.

This paper is structured as follows. In Section 2, we provide a mathematical definition of the MDCP and the MCDCP, respectively using the stochastic intensity representation. This section offers the joint Laplace transform/probability generating function, adopted from Jang and Oh (2021). The conditional expectations and variances of compound point process and their relevant conditional expectations are derived in Section 3. To do so, the martingale methodology used in Dassios and Jang (2003), based on the infinitesimal generator of the MCDCP, is employed. We also cover the derivation of joint moments of the compound point process and their linear correlations. As a measure of systemic risk, conditional value-at-risk (CoVaR) introduced by Adrian and Brunnermeier (2016) is presented in this section. As explained in Jaworski (2017), in practice, CoVaR quantifies the size of the bailout required to keep a financial institution solvent with a given probability when another financial institution incurs significant losses. For the original definition of CoVaR, we refer Adrian and Brunnermeier (2016) and its generalization in Girardi and Ergün (2013) and Mainik and Schaanning (2014).

Section 4 presents numerical computations of the moments, correlation coefficients and CoVaRs to illustrate our theoretical results on the MCDCP, where the t-copula is used for externally excited joint jumps. To compare contagious effects from catastrophic events with cases where there are no self-excited jumps, we consider the multivariate compound shot noise Cox process (MCSCP) as its counterpart. The numerical examples demonstrate that this compound point process can be used to address the interconnectedness of business and economic activities and the persistence of aftershocks across multiple business lines or the global financial system. In the context of infectious events such as cyber risk events and systemic risk events, modeling their aggregate loss with the MCDCP is well suited. We also provide the algorithm for simulating the MCDCP in this section. Section 5 concludes the paper.

2.1. Definition

Let us start with a mathematical definition for the MDCP in Definition 1 via the stochastic intensity representation. This process is within the general framework of affine process, for that see Duffie et al. (2000), Duffie et al. (2003) and Glasserman and Kim (2010). Throughout the paper, we write , ℕ₀ := ℕ ∪ {0}, ℝ₊ for the set of positive real numbers, and ℝ₀₊ for the set of non–negative real numbers. We use bold symbols for D-dimensional vectors, e.g. Nt:=(Nt(1),…,Nt(D))′ whose d^th margin is denoted with superscript (d) by Nt(d).

Definition 1. Multivariate dynamic contagion process (MDCP) is a multivariate point process N_t, t ∈ ℕ₀, where

with the non-negative−stochastic multivariate intensity process λ_t, where

with

• is a history of the joint process N_t with respect to which λ_t is adapted for t ∈ ℕ₀;

• λ₀is the vector of initial intensities at time t = 0, λ0(d)∈ℝ+;

• a is the vector of constant mean-reverting levels, a^(d) ∈ ℝ₀₊;

• δ is the vector of rates of exponential decay, δ^(d) ∈ ℝ₊ ;

• X_i, i ∈ ℕ is a sequence of i.i.d. positive externally-excited joint jumps with distribution F(x), where margins are F(x^(d)), x^(d) ∈ ℝ₊, at the corresponding random times T_1,i following a Poisson process with constant rate ρ ∈ ℝ₊, and I is the indicator function;

• Yjd(d), j_d ∈ ℕ is a sequence of i.i.d. positive self-excited jumps with distribution function G(y^(d)), y^(d) ∈ ℝ₊, at the corresponding random times T2,jd(d)generated by the intensity process λt(d);

• X_i, T_1,i, and Yjd(d)are assumed to be independent of each other.

Now we present a mathematical definition for the MCDCP in Definition 2 via the stochastic intensity representation. This compound process belongs to the more general class of affine process.

Definition 2. Multivariate compound dynamic contagion process (MCDCP) is a multivariate compound point process L_t, t ∈ ℕ₀, where

with the non-negative−stochastic multivariate intensity process λ_t which is in the form of (2.1). For d ∈ {1, 2, . . . , D},

• Ξjd(d), j_d ∈ ℕ is a sequence of i.i.d. positive individual loss amounts from risk type with distribution function J(ξ^(d)), ξ^(d) ∈ ℝ₊, at the corresponding random times T2,jd(d);

• X_i, T_1,i, Yjd(d), and Ξjd(d)are assumed to be independent of each other.

2.2. Distributional properties

As the distributional property of the MCDCP, we offer its joint Laplace transform. To do so, we start with stating the propositions adopted from Jang and Oh (2021), where L_t is MCDCP , N_t is MDCP , λ_t is multivariate shot-noise self-exciting Poisson process. Note that the propositions in Jang and Oh (2021) are based on the univariate processes, while below propositions are based on the multivariate processes.

We denote the first moment of ∏d=1DX(d) by μ_{1F^(1,2,...,D)} , the first and second moments of X^(d) by μ_1F^(d), μ_2F^(d), and the first and the second moments of Y^(d) by μ_1G^(d), μ_2G^(d), respectively. Their Laplace transforms are denoted by

where it is assumed that they are finite. We also denote the first and second moments of Ξ^(d) by μ_1J^(d), μ_2J^(d), respectively, and assume that their Laplace transforms are finite, denoted by

Proposition 1. Considering the constants, , ν_d, υ_d ∈ ℝ₀₊, and time t ∈ [0, T], the conditional joint Laplace transform and probability generating function for the multivariate intensity process λ_T , the multivariate point process N_T , and the multivariate compound point process L_T is given by

where B_d(t) is determined by non-linear ordinary differential equation (ODE)

with the boundary condition B_d(T) = υ_d and C(T) − C(t) is determined by

We can obtain the conditional joint Laplace transform for the multivariate intensity process λ_T and the multivariate compound point process L_T setting θ_d = 1 in Proposition 1.

Proposition 2. The conditional joint Laplace transform for the multivariate intensity process λ_T and the multivariate compound point process L_T given is given by

where

Setting υ_d = 0 in Proposition 2, we can obtain the conditional Laplace transform of the MCDCP L_T .

Proposition 3. The conditional Laplace transform of the MCDCP L_T given is given by

In this section, we start with deriving the conditional expectation of Lt(d), d ∈ {1, 2, . . . , D}, given initial intensity at time t = 0.

Theorem 1. The conditional expectation of the CDCP Lt(d)given λ0(d)is given by

Proof: See Appendix A.

Now we consider the conditional expectation of the product of two CDCPs, Lt(p)Lt(q) where p, q ∈ {1, 2, , . . . , D}, p ≠ q. For simplicity, without loss of generality, we choose p = 1 and q = 2 and derive the conditional expectation of Lt(1)Lt(2) directly by solving an ODE in Theorem 2, for which we start with a lemma to show the conditional expectation of λt(1)Lt(2). Note that (3.2) in Lemma 1 satisfies for any two different paths in {1, 2, , . . . , D}. Alternatively, using (2.2) it can be obtained by differentiating the Laplace transform of the MCDCP with respect to ν_d and then setting ν_d = 0. However, by solving the ODE directly it is easier to generalize to derive higher moments.

Lemma 1. The conditional expectation of λt(1)Lt(2)given λ0(1)and λ0(2)is given by

Proof: See Appendix B.

Theorem 2. The conditional expectation of Lt(1)Lt(2)given λ0(1)and λ0(2)is given by

Proof: See Appendix C.

Remark 1. Setting Af(λ,n,l,t)=Πd=1Dl(d) (A.1), we can derive E(Πd=1DLt(d)∣ℑ0). To save the space, we leave it as further research.

Based on Theorem 1 and Theorem 2, we can easily obtain the conditional linear correlation coefficient between Lt(1) and Lt(2), i.e.

We show their numerical values in Section 4. For the correlation coefficient calculation, we need the conditional variance of Lt(d), for which we start with a lemma to derive the conditional expectation of λt(d)Lt(d), d ∈ {1, 2, . . . , D}.

Lemma 2. The conditional expectation of λt(d)Lt(d)given λ0(d), is given by

Proof: See Appendix D.

Theorem 3. The conditional expectation of (Lt(d))2given λ0(d)is given by

Proof: See Appendix E.

Corollary 1. The conditional variance of CDCP Lt(d)given λ0(d)is given by

Proof: See Appendix F.

A cyber-attack may induce losses to multiple lines of business. Each external cyber risk event may generate a cluster of losses according to the branching structure of a self-exciting process. Hence the MCDCP is well suited in modeling aggregate cyber loss to capture the effect of increases of the intensity due to external cyber risk events. Its expectation and variance can be considered to calculate cyber loss insurance premiums and their numerical values are shown in Section 4.

Due to the interconnectedness of business and economic activities, exogenous shock (as occurred in the 2007–2008 global financial crisis and the COVID-19 pandemic) can trigger many following systemic failures to our financial and economic systems. Silicon Valley Bank (SVB) shut down in March 2023 is another example of showing that there exists significant risk of systemic impact to the global banking sector. Hence before closing this section, we present a generalized conditional value-at-risk (CoVaR) used in Fissler and Hoga (2023) as a measure for systemic risk and its numerical values are also shown in Section 4. The CoVaR is defined by for p₁ ∈ (0, 1) and p₂ ∈ (0, 1),

where

and we simply write CoVaRp(Lt(1)∣Lt(2))=CoVaRp∣p(Lt(1)∣Lt(2)) when p₁ = p₂ = p.

Let us provide numerical illustrations of our theoretical results on the MCDCP, i.e. the values of their moments, correlation coefficients and CoVaRs. Under three-dimensional setting, i.e. D = 3, for externally-excited joint jump distribution F(x), we use the t-copula given by Cϑ(u1,u2,u3)=tɛ,Σ(tɛ-1(u1),tɛ-1(u2),tɛ-1(u3)), where tɛ-1 is the inverse cumulative distribution function (c.d.f.) of a standard univariate t and t_ɛ,Σ denotes the c.d.f. for a multivariate t-distribution with ɛ degrees of freedom, zero mean vector, and 3 × 3 correlation matrix Σ. Here, copula parameter ϑ ∈ [−0.5, 1] is the off-diagonal element of Σ and , d ∈ {1, 2, 3}. Note that ϑ has lower bound of −0.5 to satisfy the positive definite condition of correlation matrix because we take three dimensional t-copula (see Oh et al., 2021).

For F(x^(d)), we use an exponential distribution with rate α^(d) ∈ ℝ₊. A log-gamma distribution is used for G(^(d)) and its probability density function with shape parameter c^(d) ∈ ℝ₊, rate parameter ς^(d) ∈ ℝ₊, and scale parameter ψ^(d) ∈ ℝ₊ is given by

that can capture the effect of sudden increases of intensities, i.e. aftershocks triggered by joint events.

For J(ξ^(d)), we use a generalized Pareto distribution (the beta of the second kind) whose probability density function with scale parameter ζ^(d) ∈ ℝ₊ and shape parameters ω^(d), k^(d) ∈ ℝ₊, i.e.

that can accommodate catastrophic losses due to joint events and aftershocks. We have

The log-gamma and generalized Pareto distributions are implemented using the actuar package, while the t-copula is applied using the copula package, both in R software.

To illustrate the features of the proposed process, we consider three companies that incur losses differently depending on shocks. We assume that the frequency of joint event (e.g. a contagious malware/exogenous shock) to three business lines is 3 (i.e. ρ = 3) per unit time period (say, per year) and that the contribution of the second business line to its intensity is less than that of other lines. Once the joint event occurs, it replicates itself causing a series of shockwave separately. We assume that aftershocks, which are unknown at the time of joint event’s arrival, contribute less to the intensity of the third business sector. For simplicity, the loss amount has the same distribution function for all business lines. The specific values are motivated by the Jang and Oh (2021). Hence, we set the numerical values for the parameters in Definition 1, 2 and t-copula as follows.

Table 1 presents that the expectations and variances for the MCDCP and MCSCP with D = 3 at time t = 1, 5, 10. The MCSCP has no self-excited jumps and a = (0, 0, 0). The values calculated using the MCDCP are higher than their counterparts calculated using the MCSCP, and the values for business line 1 and 3 are similar under the MCSCP. That is because three expectations of aftershocks, μ_1G^(d), d ∈ {1, 2, 3} are involved in the calculations under the MCDCP in (3.1) and (3.6). Therefore, the significance of clustering impacts of the three separate aftershocks on the intensities, driven by joint events, depends on aftershock size measures G(y^(d)), d ∈ {1, 2, 3}. This confirms that the MCDCP can be considered for modeling aggregate loss to calculate insurance premiums for contagious catastrophic events such as cyber risk events.

The conditional pairwise linear correlation coefficients between Lt(1),Lt(2), and Lt(3) at different values of ϑ in (3.3), compared to their counterparts when there are no self-excited jumps and a = (0, 0, 0), are shown in Table 2 and Figure 1. The table and figure present the correlation coefficients obtained from the equation in (3.3). For both processes, given time t, as ϑ increases, their linearities increase. We also observe that, for a given time t, Lt(1) and Lt(3) have a stronger correlation than other pairs. This is because the size of the external jumps in business lines 1 and 3 is relatively large on average, and the difference in correlation increases as ϑ increases. Moreover, the results show that the pairwise linearities computed using the MCDCP are lower than those calculated using the MCSCP at different values of ϑ. This is because the three separate aftershocks in the intensities weaken the linearities between them. Figure 2 demonstrates this, where simulations run on 20,000 paths for the MCDCP and MCSCP for t = 10, using R software and Algorithm 1.

Another interesting finding is the change of the linearity in all pairs as time t increases. As shown in Figure 1, they decrease significantly as t approaches 10 for the MCDCP, while they increase marginally for the MCSCP. These low correlations indicate that as time goes by, knowing the aggregate loss in one business line or firm does not help at all in predicting the aggregate loss in other business line or firm. Therefore, relying on correlation coefficients as a financial or economic risk measure for the MCDCP could lead to incorrect conclusions, even though the model captures the effect of intensity increases due to exogenous shocks.

Table 3 and Table 4 show CoVaRs in (3.7) as a systemic risk measure, while Table 5 shows the conditional probabilities of the aggregate loss in the first business line exceeding its VaR, given that the aggregate loss in the second business line exceeds its VaR,

The values are obtained under various scenarios, which are combination of t = 1, 5, 10, ϑ = −0.5, 0, 0.5, 0.99, and . Note that the values in Tables 3–5 are mean values over 50 iterations, with 20,000 paths simulated for each iteration. For both the MCDCP and the MCSCP, given t and , the risk measures tend to increase as ϑ increases. We observe that for each t, ϑ, , the values of CoVaRs calculated using the MCDCP are higher than those calculated using the MCSCP. In contrast, the conditional probabilities calculated using the MCSCP are higher than those calculated using the MCDCP , consistent with the correlation values in Table 2.

Both risk measures, CoVaRs and conditional probabilities, focus on the tail of the joint distribution but exhibit opposite features. This confirms the appropriateness of CoVaR as a systemic risk measure for the MCDCP when dealing with tail risks for infectious events. Additionally, employing a broader mix of risk measures, rather than relying on a single risk measure, could help avoid the pitfalls of making drastically erroneous decisions, especially when weak linearities exist between the risks involved in practice.

For a malicious virus or exogenous shock, it is critical the spread be isolated from the surrounding network, server, or financial system as quickly as possible, lest it rapidly multiply in scale. The U.S. government stepped in to guarantee deposit accounts at SVB in March 2023 to restore confidence in the global banking sector. Public and private sector organizations are focused on enhancing IT and network security by implementing effective cybersecurity plans and strategies in place. The moment values of the MCDCP in Table 1 highlight that the aggregate loss increases over time, and this rise is much more significant than conceptualized under previous models that do not account for the contagious aftershocks beyond an initial outbreak. Properly accounting for this effect allows for a more accurate understanding of the risks posed by of malicious viruses or exogenous shocks, with potential applications in pricing a market which can adequately insure these risks and in designing more effective central bank or government intervention strategies.

We close this section providing the simulation algorithm for one sample path of the MCDCP (L_t, N_t, λ_t) truncated at time T to make it easier for statistical analysis, further business applications and research (see Figure 3). For the simulation of the univariate CDCP, we refer Jang and Oh (2021) Section 4 algorithms.

We have introduced the MCDCP and presented its analytical expression for the joint Laplace transform and probability generating function. Under the MCDCP, we have studied the correlation coefficients, moments and conditional value-at-risks (CoVaRs). These moments can be used to compute actuarial net and gross premiums.

The proposed MCDCP is designed to address the interconnectedness of business and economic activities and the persistence of aftershocks between the risks involved and their contagious loss aggregation. It can serve as a model for enterprises, organizations, central banks, and governments to quantify contagious catastrophic losses across multiple business lines or the global financial system. Modeling aggregate loss from infectious events with the MCDCP is well-suited for contexts such as cyber risk events and systemic risk events. Since the proposed MCDCP accounts for simultaneously occurring joint events across multiple business lines or the global financial system via multivariate intensity process, it can be extended in the future by incorporating dependence structure modeling in the infinitesimal generator of the MCDCP, which we leave as the next objective for further research. We have also provided an algorithm for simulating the MCDCP, which can be used for statistical analysis, further business applications, and research purposes.

Given the broad applicability of the MCDCP, we anticipate that what we have presented in this paper will provide practitioners with feasible approaches to quantify multivariate risks across a range of domains, including disruptive technologies, environmental changes, global pandemics, ecology and epidemiology, as well as in the business and economic sectors. The next step in fitting the MCDCP involves obtaining the event arrival time points for spillover events in these fields, which we also leave for further research.