Cyber risk modeling: a discrete multivariate count process approach

Abstract

In the past decade, cyber risk has raised much interest in the economy, and cyber risk has evolved from a type of pure operational risk to both operational and liability risk. However, the modeling of cyber risk is still in its infancy. Compared with other financial risks, cyber risk has some unique features. In particular, discrete variables regularly arise both in the frequency component (e.g. number of events per unit time), and the severity component (e.g. the number of data breaches for each cyber event). In addition, the modeling of these count variables are further complicated by nonstandard properties such as zero inflation, serial and cross-sectional correlations, as well as heavy tails. Previous cyber risk models have largely focused on continuous models that are incompatible with many of these characteristics. This paper introduces a new count-based frequency-severity framework to model cyber risk, with a dynamic multivariate negative binomial autoregressive process for the frequency component, and the generalized Poisson inverse-Gaussian distribution for the severity component. We unify these new modeling tools by proposing a tractable Generalized Method of Moments for their estimation and applying them to the Privacy Rights Clearinghouse (PRC) dataset.

Keywords:

• Operational risk
• cyber risk
• systemic and pairwise contagion
• self-excitation
• negative binomial autoregressive process
• generalized poisson inverse-Gaussian distribution
• discrete stable distribution

JEL classifications:

• C22
• G17

Acknowledgments

The authors acknowledge helpful comments and suggestions from Editor Mogens Steffensen and an anonymous referee.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1 Currently, there is no consensual definition of cyber risk in the insurance sector. However, a few associations of insurance companies have put forward some definitions. For example, the Geneva Association (2016) defined cyber risk as ‘any risk emerging from the use of information and communication technology that compromises the confidentiality, availability, or integrity of data or services.’ Similarly, Cebula & Young (2010) define cyber risk as ‘operational risks to information and technology assets that have consequences affecting the confidentiality, availability or integrity of information or information systems’.

2 Traditionally, operational risk is assumed to be either uncorrelated, or negative correlated with other risks (Gatzert & Kolb 2014).

3 A recent literature introduces network-based (Fahrenwaldt et al. 2018, Jevtić & Lanchier 2020) or epidemiology-based (Xu & Hua 2019, Hillairet & Lopez 2021) frequency models. Even though these models have appealing theoretical properties, they are not suitable for currently available cyber risk databases. As Hillairet & Lopez (2021) put it, ‘Using the model to replicate the attack is a difficult task, since few available data support the calibration of the parameters.’ As a consequence, the illustration of these models have been largely confined to simulation studies.

4 The Privacy Rights Clearinghouse database is made publicly available at https://privacyrights.org/data-breaches. This article is based on a download on 2/27/2021, covering events between 2005 and 2019.

5 Eling & Jung (2018) also suggest a categorization according to organization types. Our empirical analysis in Section 5 will focus on classification according to breach types due to length restriction. We find consistent results across different organization types, and we collect the estimation results for organization types in Appendix 3.

6 The weekly frequency has also been adopted by Jung (2021). By using weekly counts, we obtain a much longer time series, which hugely facilitates the parameter estimation of our model. Similarly, higher-frequency count data have also been analyzed in the cyber risk literature by, e.g. Xu et al. (2017) and Peng et al. (2017). Other choices of frequencies include bi-weekly (Jung 2021), monthly (Wheatley et al. 2016, Eling & Jung 2018, Eling & Wirfs 2019, Xu & Zhang 2021), bi-annually (Wheatley et al. 2021), and yearly (Wheatley et al. 2016).

7 Depending on the applications, this unit might be one year, in the case of computation of Value-at-Risk for a bank (Chavez-Demoulin et al. 2005, 2016), or one hour (Peng et al. 2017).

8 Bessy-Roland et al. (2021) show that the computation of the unconditional expectation of the process involves matrix equations of dimension $K(K+1)\times K(K+1)$, where K is the dimension of the Hawkes process. The conditional expectation given observation up to a certain time t has not been derived.

9 Some of the commercial databases, such as SAS OpRisk Global data (e.g. Eling & Wirfs 2019, Eling & Schnell 2020), or IBM's annual Cost of Data Breach Report (Eling & Jung 2018), also provide the monetary cost of these cyber events. However, these datasets are not publicly available and most existing cyber risk studies rely on data that provide integer-valued severity variables. See also Section 6 and Appendix 2 for further discussion.

10 According to the RMS-CCRS review, 46% of the reviewed policies included coverage for reputation damage (RMS-CCRS 2016).

11 For instance, on 9/13/2016, the World Anti-Doping Agency (WADA) confirmed the hack of its database on medical files of Olympic athletes by the Russian group Fancy Bear. This event led to only 4 lost records, but its economic consequence could be difficult to evaluate, yet potentially large, especially if the stolen information was used by rival athletes.

12 Trend in cyber data has been discussed in the literature. For example, Edwards et al. (2016) report no time trend for the severity over the past decade. Eling & Loperfido (2017) report no clear upward or downward trend in the PRC data from 10 January 2005 to 15 December 2015. According to Eling & Jung (2018), they empirically do not observe any non-stationarity or serial dependence for different risk types. Jung (2021) reports a break in loss severity trend using 50 days rolling windows.

13 For instance, if we were to count the cyber events in PRC by organization types as in Eling & Jung (2018), then we get K = 5. These results are collected in Appendix 3.

14 Under some conditions, the domain of ${\delta }_0$ can be further relaxed, see Denuit & Lu (2021) for a discussion.

15 Even though the covariance formula (2(2) $\mathrm{Corr}[X_i,X_j]=\frac{{\sigma }_{\mathrm{ij}}^2}{{\sigma }_{\mathrm{ii}}{\sigma }_{\mathrm{jj}}}\ge 0,\mathrm{\forall }i,j=1,\dots ,K.$(2) ) does not depend on the sign of the entries Σ, the off-diagonal elements of Σ should not be constrained to be nonnegative. Indeed, the distribution of $X$ is not multivariate normal, hence it is not fully-characterized by its first two moments.

16 This latter involves quadratic forms of bivariate Wishart stochastic matrices. The problem with this specification is that at dimension K, each such quadratic form involves $K(K+1)/2$ parameters, and there are as many quadratic forms as there are components of $Y_t$. In contrast, our new model only involves the Wishart-gamma distribution, which is the diagonal of the Wishart distribution and is hence much more parsimonious.

17 The effect of self-excitation is captured by parameters ${\beta }_{j,j},j=1,\dots ,K$.

18 The effect of mutual-excitation is captured by parameters ${\beta }_{i,j},i\ne j.$

19 Note that the term ‘systemic contagion’ does not necessarily mean that one systemic event generates simultaneously cyber events of different types. This terminology is mainly introduced to emphasize the fact that in contrast to self-excitation which generates counts that are conditionally independent given the past, the systemic channel, through the common factor $X_t$, captures the positive correlation of $Y_{1,t+1},\dots ,Y_{K,t+1}$ given past information.

20 See Chavez-Demoulin et al. (2016) for a similar, deterministic approach for the modeling of operational risk.

21 We say that process $(Y_t)$ is strictly stationary, if any finite dimensional marginal distribution of $(Y_{t+k_1},Y_{t+k_2},\dots ,Y_{t+k_n})$, where $n,k_1,\dots ,k_n$ are integers, are time-invariant.

22 See also El-Shaarawi et al. (2011) for an application in environmental sciences.

23 See Zhu & Joe (2009) for the interpretation of these three parameters in terms of shape, scale, and probability, respectively.

24 Recently, Sun et al. (2020) investigate the dependence between the severity and frequency at the firm-level and found a nonlinear, positive dependence. However, their analysis is based on firms that have reported claims only, and is hence subject to selection bias. Moreover, their approach does not allow for serial correlation, and cannot be applied at the portfolio level. Because events are disaggregated into different breach categories, some of the heterogeneity has already been taken into account. As a consequence, the independence assumption between frequency and severity for each of the event types seems reasonable.

25 Zhu & Joe (2009) propose an alternative recursive formula to compute the pmf, but the computational cost of this recursion is still very high for large n and it is known that such recursion may have large numerical error (Wang & Panjer 1993).

26 Note, however, that it is possible to compute the conditional pmf of a multivariate count process by series expansion, if the conditional pgf is available in closed form. See Lu (2021) for a treatment of the univariate case, and Darolles et al. (2019) for an example in dimension K = 2. The disadvantage of this approach is that the computation becomes much more costly in high dimensions, say K = 4, even though one could rely on pairwise conditional pgf $p(Y_{t,i},Y_{t,j}|Y_{t-1})$, $1\le i\ne j\le 4$, and estimate the model using (pairwise) composite likelihood approach (see Pedeli & Karlis 2013, for an example in multivariate count time series analysis). As a consequence, in this paper, we focus on the GMM approach, which is much more computationally tractable.

27 Then the corresponding standard errors of the estimator $\hat{\mathrm{\Sigma }}$ are calculated through the Delta Method.

28 Note, that even though the log-likelihood function has a closed form, in the estimation section, we did not use this result for estimation, but resorted to GMM instead. This is due to the fact that the GMM is much more computationally feasible for the estimation of our multivariate NBAR process.

29 Indeed, it suffices to adapt the proof of Proposition 3, by taking second order derivative in Equation (A3(A3) $\mathbb{E}[W^{1-\beta }\exp (-\mathrm{tW})]=\frac{1}{\mathrm{\Gamma }(\beta )}{\int }_0^{\mathrm{\infty }}{s}^{\beta -1}\mathbb{E}[e^{-(t+s)W}W]\mathrm{d}s,$(A3) ) instead of the first order derivative.