# Incident-Specific Cyber Insurance

Wing Fung Chong  Daniël Linders  Zhiyu Quan  Linfeng Zhang Corresponding
author.

###### Abstract

In the current market practice, many cyber insurance products offer a coverage
bundle for losses arising from various types of incidents, such as data
breaches and ransomware attacks, and the coverage for each incident type comes
with a separate limit and deductible. Although this gives prospective cyber
insurance buyers more flexibility in customizing the coverage and better
manages the risk exposures of sellers, it complicates the decision-making
process in determining the optimal amount of risks to retain and transfer for
both parties. This paper aims to build an economic foundation for these
incident-specific cyber insurance products with a focus on how incident-
specific indemnities should be designed for achieving Pareto optimality for
both the insurance seller and buyer. Real data on cyber incidents is used to
illustrate the feasibility of this approach. Several implementation
improvement methods for practicality are also discussed.

Keywords: Risk management, Cyber insurance, Incident specificity, Statistical
learning, Pareto optimality.

##  1 Introduction

In May 2021, Colonial Pipeline, one of the largest oil pipeline systems in the
United States (US), was hit by a ransomware attack and had a shutdown for five
days; see [22]. In November 2020, Amazon Web Services, a major cloud service
provider that many businesses rely on, had a severe outage triggered by an
operating system configuration; see [13]. In 2017, Equifax, one of the three
largest US credit reporting agencies, experienced a major data breach that
exposed the private records of nearly 150 million American citizens; see [10].
In 2016, a class action was brought against Meta Platforms (formerly known as
Facebook) for failing to comply with the Biometric Information Privacy Act in
Illinois, and the case was settled for $650 million in 2021; see [14]. These
examples illustrate the extraordinary exposure of modern businesses to cyber
risk and their possible devastating consequences.

The events listed above represent a wide range of incidents that are
considered cyber-related, and they are common in this digital world where
business operations rely heavily on data and cyber systems. To mitigate the
aftermath impact of such incidents, cyber insurance products are designed for
businesses to purchase. The insurability of cyber losses is examined in [8],
in which their loss expectations are shown to be finite. These cyber insurance
products typically cover losses resulting from various perils. As summarized
in [23, 19, 27], some cyber insurance policies indemnify losses caused by
ransomware, and some policies cover liabilities resulting from data breaches.
It is also not unusual for a single policy to cover multiple types of cyber
incidents. In this paper, our focus is on those policies that offer coverage
for multiple perils, which are referred to as incident-specific cyber
insurance.

### Need for incident-specific cyber insurance

Incident-specific cyber insurance products are in line with businesses’ needs
for cyber risk management.

A definition of cyber risk commonly adopted in practice and literature is
“cyber risks are operational risks to information and technology assets that
have consequences affecting the confidentiality, availability, and integrity
of information and information systems”; see [6]. Aligning to this idea, the
National Institute of Standards and Technology (NIST) Cybersecurity Framework,
from which many cybersecurity standards for critical infrastructure sectors
are derived, considers data security from the perspectives of data
confidentiality, integrity, and availability and proposes a comprehensive set
of best practices; see [21]. In [9], the authors argued that cyber risks are
subcategories of operational risk and used cyber-related data from a database
of operational losses to model the loss of four types of cyber incidents,
including the ones caused by actions of people, system failure, internal
process failure, and external events. In [2], the author considered the same
four categories of cyber risks and used Bayesian networks to model the
relationship between risk categories and risk factors. In [11], the authors
proposed a different set of incident types for supply chain cyber risk
management, including physical threats, breakdown, indirect attacks, direct
attacks, and insider threats.

Provided that cyber risk is a collective term referring to the risks of
multiple types of cyber incidents, and businesses typically demand a risk
management strategy that is holistic and, at the same time, can address the
subtle differences among various incident types, purchasing an incident-
specific cyber insurance policy is then a natural choice for businesses to
meet such a need.

### Existing incident-specific cyber insurance

To ensure risk exposures are manageable, many incident-specific cyber
insurance policies come with sublimits and separate deductibles for individual
coverages. For example, [1] allows for setting limits and deductibles
individually for events such as cyber extortion and computer crime. [23]
surveyed 67 policies and found many of them use sublimits. In this paper,
these sublimits and deductibles shall be referred to as incident-specific
limits, or simply limits, and incident-specific deductibles or deductibles.
The scope of this study though does not include the cases in which there is a
limit or deductible at the aggregate policy level, which shall be deferred to
future research.

Although insurers can shield them from excessive exposures using those risk-
sharing provisions, such policy designs complicate the process of shopping for
cyber insurance coverage for prospective buyers. Given that most companies
have difficulties in determining the appropriate amount of coverage (see
[15]), this problem will only be more of an obstacle that hinders cyber
insurance purchases if companies have to impromptu configure the amounts of
coverage for individual incident types.

### Proposed workflow for determining coverage

Provided that incident-specific cyber insurance is a norm in practice but
lacks a technical foundation that justifies the determination of coverage
amounts, it is in the interests of both insurers and companies, who seek
insurance coverage, to make the process of determining incident-specific
coverage easy to understand, compatible with the existing underwriting
procedures, and mutually beneficial to both parties. To this end, we shall
address the aforementioned problem of determining incident-specific coverage
from a bilateral perspective.

Drawing inspiration from [4], which shows that an environment-specific
indemnity profile could be Pareto optimal for both the buyer and seller, in
this paper, we propose an economically sound workflow that determines the
appropriate amount of incident-specific coverage. By doing so, it helps the
insured and the insurer reach a Pareto optimality, such that the risk taken by
either of the two parties cannot be further reduced without increasing the
other’s risk. How such a problem should be formulated and how both parties’
risks are measured are elaborated in Section 2.

The diagram in Figure 1 illustrates this workflow, which consists of three
major parts: (I) the incident-specific cyber insurance problem, which serves
as the central pillar of this workflow and assigns optimal amounts of risk to
the policyholder and the insurer; (II) model inputs, in which severity and
incident type can be estimated and statistically learned based on historical
cyber incident data and the underwriting process, together with the
policyholder’s and insurer’s preferences information for the decision-making
part; and (III) the solver for the optimal insurance problem, which can be
implemented using either an exact method or an approximated but quick method
depending on whether there is a time constraint. The solver produces
specifications of the incident-specific insurance coverage, which is the final
and desired outcome of this workflow.

Figure 1: Workflow of designing incident-specific cyber insurance coverage

Regarding Part (III), a practical scenario is that the insurer needs to
provide a quote and a coverage to a prospective customer on-the-fly. We shall
show that using the exact method may result in a long waiting time, and thus
the quick method, which generates the results almost instantly, can be a
workaround in that situation.

In this paper, we shall detail the key components of the proposed workflow for
designing incident-specific cyber insurance coverage, and then demonstrate it
with numerical examples.

### Findings and contributions

First, based on the incident-specific insurance design problem, the cyber
insurance contract supplied by the proposed workflow must be reducing the
aggregate cyber risk of a company, as an insurance buyer, and the insurer
after the risk transfer; the contract must also be mutually beneficial to both
the insured and the insurer since a paid premium is less than the buyer’s
cyber risk reduction while is able to cover the risk taken by the insurer.
However, their aggregate cyber risk reduction is not necessarily effective, in
the sense that some insurance buyers might pay more premium than their
reduction in the aggregate risk with the insurer; this is indeed observed from
the numerical examples presented in Section 5. Second, the proposed workflow
is computationally tractable in scenarios that require the specifications of
an incident-specific policy to be delivered promptly, such as quick quote and
coverage services, comparing to other conventional methods.

Therefore, the contributions of this paper are as follows. Firstly, cyber
insurance policies commonly cover multiple perils, but the determination of
incident-specific coverage amounts lacks justification in practice and
literature. The proposed workflow leads to insurance specifications that are
justified and economically sound to both the insured and the insurer.
Secondly, the proposed workflow is data-driven and compatible with the
existing underwriting procedures. All the data required by this workflow can
be obtained from public or proprietary sources by underwriters, which results
in no additional work on the frontend. Thirdly, we show that, from the
numerical examples, while the designed incident-specific cyber insurance
contract must be mutually beneficial to the insured and the insurer, their
aggregate cyber risk reduction might not be always effective. Lastly, this
paper overcomes the possible computational challenges arising from the
implementation of such a workflow, making it versatile and adaptable to
different scenarios depending on the amount of time permitted to produce
results.

The rest of this paper is organized as follows. Section 2 formulates the
problem of designing incident-specific cyber insurance as an optimization
problem. Section 3 discusses the method and the data that we used to obtain
the necessary model inputs. Given the model inputs, Section 4 describes the
Cross Entropy Method and how it is used to solve the proposed problem as well
as presents and discusses some numerical results. Section 5 presents the use
of function approximation to address some computational challenges in the
process of solving the optimization problem. We conclude in the last section
with discussions of potential applications and future works.

##  2 Design of Incident-Specific Cyber Insurance

###  2.1 Buyer and seller’s Pareto optimality

Incident-specificity in this study means that although the insurance product
is designed as a package with coverage for a variety of cyber incident types,
the coverage for each incident type has its own indemnity function. In
addition, a major characteristic of the policy design is that different types
of cyber incidents are mutually exclusive, i.e., every incident precisely
belongs to a unique category. Therefore, the occurrence of one incident would
only trigger one and only one corresponding incident-specific coverage. This
point will be further illustrated with real data in later sections. Then, the
problem of interest in this study is how to determine the incident-specific
indemnity functions.

To formulate this problem, let $(\Omega,\mathbb{F})$ be a measurable space,
let $\mathbb{P}$ be a probability measure on $(\Omega,\mathbb{F})$, and let
$\mathcal{K}=\left\\{1,2,\dots,K\right\\}$ be the set of indices of insurable
incident types with cardinality $K$. For a cyber incident of any type, denote
the random variable of the realized non-negative loss as
$(X)\mathbb{I}_{\\{O=k\\}}$, where $X$ is the loss resulting from a particular
cyber event and the indicator function $\mathbb{I}_{\\{O=k\\}}$ takes the
value of $1$ if the Occurred incident type $O$ is $k$, for any
$k\in\mathcal{K}$, and takes the value of $0$ otherwise. Furthermore, let
$I_{k}$ be the incident-specific Indemnity function, for incident type
$k\in\mathcal{K}$; that is, the insurance buyer receives
$I_{k}(X)\mathbb{I}_{\\{O=k\\}}$ after an incident. The amount of loss
Retained by the buyer is
$R_{k}(X)\mathbb{I}_{\\{O=k\\}}=(X)\mathbb{I}_{\\{O=k\\}}-I_{k}(X)\mathbb{I}_{\\{O=k\\}}$,
which could be non-zero in the scenario that an incident-specific indemnity
does not fully cover the corresponding loss.

Acquiring such an incident-specific insurance policy costs the buyer a single
premium $\pi$. From the buyer’s perspective, the total loss at the end of the
policy period is, $\pi+\sum_{k=1}^{K}R_{k}(X)\mathbb{I}_{\\{O=k\\}}$, where
the summation is a result of the mutual exclusivity of incident types. The
risk-sharing counterparty, i.e., the seller, correspondingly bears the risk of
value $\sum_{k=1}^{K}I_{k}(X)\mathbb{I}_{\\{O=k\\}}$ in exchange for receiving
the premium $\pi$, thus making the loss at the end of the policy period,
$-\pi+\sum_{k=1}^{K}I_{k}(X)\mathbb{I}_{\\{O=k\\}}$. In this study, any
potential return generated by the premium during the policy period is
neglected for simplicity.

To support the decision-making regarding configurations of the contract
consisting of the incident-specific indemnity functions and the policy
premium, assume that the Buyer and the Seller are endowed with risk measures
$\rho_{B}$ and $\rho_{S}$, respectively, and the Pareto optimality is to be
achieved with respect to the risk measures on their post-transfer risk
positions, so that neither of the two parties can attain a lower risk,
evaluated by the risk measure, without the other party’s risk being increased.
Then, as discussed in [4] and proved in Theorem 3.1 of [3], if their risk
measures are translation invariant (which are satisfied by most risk
measures), their Pareto optimality is achieved by the contract which
minimizes:

$F(I_{1},I_{2},\dots,I_{K},\pi;X,O)=\rho_{B}\left(\sum_{k=1}^{K}R_{k}(X)\mathbb{I}_{\\{O=k\\}}+\pi\right)+\rho_{S}\left(\sum_{k=1}^{K}I_{k}(X)\mathbb{I}_{\\{O=k\\}}-\pi\right).$

With their risk measures being translation invariant and the policy premium
$\pi$ being a deterministic variable, the objective function, with a slight
abuse of notation, can be simplified to

$F(I_{1},I_{2},\dots,I_{K};X,O)=\rho_{B}\left(\sum_{k=1}^{K}R_{k}(X)\mathbb{I}_{\\{O=k\\}}\right)+\rho_{S}\left(\sum_{k=1}^{K}I_{k}(X)\mathbb{I}_{\\{O=k\\}}\right),$

The minimization problem that solves the ex-ante specified Pareto optimal
incident-specific indemnity functions is then given by:

$\min_{(I_{1},I_{2},\dots,I_{K})\in\mathcal{I}}F(I_{1},I_{2},\dots,I_{K};X,O),$
(1)

where, with $\mathrm{Id}$ being the identity function, the solution space

$\mathcal{I}:=\\{(I_{1},I_{2},\dots,I_{K}):0\leq
I_{k}\leq\mathrm{Id},\;I_{k}\text{ and }R_{k}\text{ are non-
decreasing},\;\forall k\in\mathcal{K}\\}$

adheres to two fundamental principles in insurance, including 1) the buyer is
compensated partially or fully for a loss, but cannot make a profit from the
compensation; and 2) the ex-post moral hazard that a falsely larger claim is
made to reduce the buyer’s deductible or the seller’s indemnity payout should
be avoided.

Finally, the ex-ante specified Pareto optimal policy premium is given by the
rationality constraints of the buyer and seller:

$\rho_{B}\left(\sum_{k=1}^{K}R_{k}(X)\mathbb{I}_{\\{O=k\\}}+\pi\right)\leq\rho_{B}\left(X\right),$
$\rho_{S}\left(\sum_{k=1}^{K}I_{k}(X)\mathbb{I}_{\\{O=k\\}}-\pi\right)\leq\rho_{S}\left(0\right)=0,$

in which the right-hand-side of the inequalities are the risk measures on both
parties’ pre-transfer risk positions. By the translational invariance of their
risk measures, these can be simplified as:

$\rho_{S}\left(\sum_{k=1}^{K}I_{k}(X)\mathbb{I}_{\\{O=k\\}}\right)\leq\pi\leq\rho_{B}\left(X\right)-\rho_{B}\left(\sum_{k=1}^{K}R_{k}(X)\mathbb{I}_{\\{O=k\\}}\right).$
(2)

in which the incident-specific indemnity and retained functions are Pareto
optimal being solved in (1). Thus, at the Pareto optimality, the policy
premium is not disentangled from the respective risks taken by the buyer and
the seller. Varying the Pareto optimal policy premium choice within this
interval traces the entire Pareto frontier. Note that, in this paper, the
premium is derived from an economic perspective; in the insurance practice, an
insurer could charge a premium based on an actuarial price.

Therefore, in the remainder of this study, the focus is on the design of the
incident-specific indemnity functions.

###  2.2 Optimal indemnities with Value-at-Risk preferences

The choice of risk measures, $\rho_{B}$ and $\rho_{S}$, depends on the
insurance buyer’s and seller’s own risk management appetites. Some common
choices, which satisfy the translational invariance, include Value-at-Risk
(VaR) and Tail Value-at-Risk (TVaR).

It can be verified that if both $\rho_{B}$ and $\rho_{S}$ are the same type of
subadditive risk measures and only differ in their parameters (e.g., TVaR with
different risk tolerance level parameters), then the solution to Problem (1)
assigns all risks to the party that has a higher risk tolerance. VaR, on the
other hand, is not subadditive but is still prevalently adapted in the banking
and insurance industries; see, for example, the capital requirement in Basel
III and Solvency II frameworks. Therefore, in this study, we assume that the
buyer and the seller both adopt VaR as their risk measures. VaR is as follows:

$\text{VaR}_{\gamma}(Y)=\inf\\{y\in\mathbb{R}:\mathbb{P}(Y\leq
y)\geq\gamma\\},$ (3)

where $Y$ is a random variable and $\gamma\in(0,1)$ is the risk tolerance
level.

With the VaR risk preferences, [4] showed that the indemnity functions in the
sub-solution space $\mathcal{I}\setminus\mathcal{I}_{1}$ are at least
suboptimal for Problem (1), where

$\displaystyle\mathcal{I}_{1}:=\\{(I_{1},I_{2},\dots,I_{K})\in\mathcal{I}:$
$\displaystyle\;I_{k}(X)=(X-d_{k})_{+}\text{ or }I_{k}(X)=X-(X-d_{k})_{+},$
(4) $\displaystyle\;d_{k}\in[0,\operatorname*{ess\,sup}(X)],\text{ for each
}k\in\mathcal{K}\\},$

where $\operatorname*{ess\,sup}(X)$ is the essential supremum of $X$ under the
probability measure $\mathbb{P}$. This implies that, for each incident-
specific coverage, either a deductible or a policy limit, denoted by $d_{k}$,
should be implemented. Within the solution space in (4), solving the infinite-
dimensional Problem (1) is thus reduced to solve the following finite-
dimensional, but combinatorial, minimization problem:

$\min_{\begin{subarray}{c}\bm{d}\in\mathbb{R}_{+}^{K},\\\
\bm{\theta}\in\\{0,1\\}^{K}\end{subarray}}F(\bm{d},\bm{\theta};X,O)=\min_{\begin{subarray}{c}\bm{d}\in\mathbb{R}_{+}^{K},\\\
\bm{\theta}\in\\{0,1\\}^{K}\end{subarray}}\left[\text{VaR}_{\alpha}\left(L_{S}(\bm{d},\bm{\theta};X,O)\right)+\text{VaR}_{\beta}\left(L_{B}(\bm{d},\bm{\theta};X,O)\right)\right],$
(5)

where

  * •

$\bm{\theta}=(\theta_{1},\theta_{2},\dots,\theta_{K})$, and for each
$k\in\mathcal{K}$, $\theta_{k}\in\left\\{0,1\right\\}$ is the choice between a
limit or a deductible being implemented. Here, we set $\theta_{k}=0$ if an
incident-specific limit is placed, and $\theta_{k}=1$ if a deductible is
implemented instead;

  * •

$\bm{d}=(d_{1},d_{2},\dots,d_{K})$ are the amounts of incident-specific limits
or deductibles, depending on which of the two is imposed on the coverage for
each incident type;

  * •

$\alpha$ and $\beta$ are choices of risk tolerance levels of the seller and
the buyer, respectively;

  * •

$L_{S}(\bm{d},\bm{\theta};X,O)=\sum_{k=1}^{K}\Big{(}\theta_{k}\left(X-d_{k}\right)_{+}+(1-\theta_{k})\left(X-\left(X-d_{k}\right)_{+}\right)\Big{)}\mathbb{I}_{\\{O=k\\}},$
(6)

and

${L_{B}(\bm{d},\bm{\theta};X,O)=\sum_{k=1}^{K}\Big{(}(1-\theta_{k})\left(X-d_{k}\right)_{+}+\theta_{k}\left(X-\left(X-d_{k}\right)_{+}\right)\Big{)}\mathbb{I}_{\\{O=k\\}}},$
(7)

are the seller’s and the buyer’s loss random variables, respectively, after
agreeing with an insurance contract.

To solve Problem (5) for Pareto optimal $\bm{\theta}$ and $\bm{d}$, the
probability distributions of $L_{S}$ and $L_{B}$ in (6) and (7) are necessary
model inputs to calculate the VaRs in the objective function; these are
discussed in Section 3. Though Problem (5) is of finite-dimensional in
$\bm{d}$, it is also of combinatorial type such that the number of
combinations in $\bm{\theta}$ grows exponentially in $K$, making the
minimization problem not mathematically tractable; see the number of cases to
consider in Proposition 3.1 of [4] which seeks for explicit optimal
indemnities even if $K=2$. Section 4 discusses a numerical method and its
algorithm to solve the combinatorial optimization problem.

##  3 Model Inputs

The distributions of the seller’s and the buyer’s loss random variables in (6)
and (7) belong to the mixture class. Therefore, model inputs to solve Problem
(5) are (i) the index set $\mathcal{K}$ of insurable cyber incident types,
(ii) the probabilities, $\bm{p}=(p_{1},p_{2},\dots,p_{K})$, of an incident
being in an individual category, where $p_{k}={\mathbb{P}(O=k)}$, for
$k\in\mathcal{K}$ (note that, due to the mutual exclusivity of incident types,
$\sum_{k=1}^{K}p_{k}=1$), and (iii) the incident-specific severity of the
ground-up loss $(X)\mathbb{I}_{\\{O=k\\}}$, for $k\in\mathcal{K}$.

A dataset from Advisen consisting of historical cyber incidents is used to
obtain the inputs to determine the loss random variables. In the dataset,
there are 103,061 historical cyber incidents. Its earliest observation dates
back to 1987, and the latest observation was recorded in September 2018. This
dataset has 125 explanatory variables coded for each incident, which can
generally be categorized into four groups as follows: (i) nature of the
incident, (ii) victim company information, (iii) consequences of the incident,
and (iv) information from any associated lawsuits. With the information on
incident types and loss amounts, we can build statistical models to predict
the probability vector $\bm{p}$ and the conditional severity distribution of
$(X)\mathbb{I}_{\\{O=k\\}}$, for $k\in\mathcal{K}$. The data and models shall
be further elaborated on in the following subsections.

###  3.1 Types of cyber incident

This paper focuses on the cyber incidents that took place in the US, which
represent close to 80% of all observations. Based on the country information
provided for each incident in the dataset, we extract a sample of all US-based
incidents, and the sample size is 84,938.

The response variable is the cyber incident type. The different incident types
are the following: (1) cyber extortion, (2) data - malicious breach, (3) data
- physically lost or stolen, (4) data - unintentional disclosure, (5) denial
of service (DDOS)/system disruption, (6) digital breach/identity thief, (7)
identity - fraudulent use/account access, (8) industry controls & operations,
(9) IT - configuration/implementation errors, (10) IT - processing errors,
(11) network/website disruption, (12) phishing, spoofing, social engineering,
(13) privacy - unauthorized contact or disclosure, (14) privacy - unauthorized
data collection, and (15) skimming, physical tampering. Similar to [17], in
this paper, these 15 categories are grouped into four types, which are (1)
data breach, (2) fraud and extortion, (3) IT error, and (4) privacy violation.
The following provides a brief description of each type of cyber incident
classified and, together with the “other” level, Table 1 shows the number of
observations for each type of cyber incident classified. Abbreviations in
parentheses will occasionally be adopted throughout this paper for clear
exposition.

  * •

Privacy violation (PV) incidents occur when companies collect or disclose
individuals’ sensitive information, such as personally identifiable
information and financial information, without receiving consent from those
individuals.

  * •

Data breaches (DB) are incidents in which data storage devices are breached,
causing a possible leakage of confidential information. These incidents can be
caused by either hacking activities or the loss of physical devices.

  * •

Cyber incidents in the fraud and extortion (FE) category are similar to
traditional frauds and extortion events, but take place in cyberspace. Cyber
frauds typically have forged digital identity involved, such as phishing
attacks, and commonly in cyber extortion events, information and information
systems are held hostage by intruders for financial gain.

  * •

In IT error (ITE) events, there is no malicious intent involved. They are
caused by incorrectly configured or operating IT systems.

Incident Categories | Number of Observations
---|---
(PV) Privacy Violation | 51315
(DB) Data Breach | 26492
(FE) Fraud and Extortion | 4464
(ITE) IT Error | 2102
Other | 565
Table 1: Number of historical cyber incidents in each category.

###  3.2 Incident type occurrence probabilities

Company-specific information is useful when estimating the occurrence
probabilities of different incident types. For example, it is reasonable to
assume that technology companies with lots of customer data are more likely to
experience data breaches than manufacturing companies, whereas manufacturers
with complex industrial control systems are more prone to IT errors that cause
disruptions than technology companies. Therefore, we use company
characteristics such as industry as predictors for the occurrence
probabilities. In this subsection, we present several predictive models that
we experimented with for this purpose and make comparisons among them. This
approach is similar to insurance rate-making which determines insurance
premiums based on the risk characteristics of the insurance buyers.

#### Selected explanatory variables

Among all explanatory variables available in this dataset, we desire
manageable and meaningful ones, and thus the variable selection process
incorporates the following considerations.

  * •

The values of some variables can only be observed after the occurrence of an
incident, at which point the type of incident is already known. For example,
the settlement cost of a lawsuit associated with an incident is not observable
when the incident occurs, and thus this cost has no impact on the type of the
incident. Such variables are ruled out from the models.

  * •

Categorical explanatory variables which have a large number of levels and are
downstream in the hierarchy, if any, are removed. For example, among the
explanatory variables that represent the geographical information of a
company, the one that contains city information is removed, but the one with
state information is kept. Likewise, explanatory variables that are more
granular industry classification codes, such as Standard Industrial
Classification (SIC) codes with more than 2 digits, are removed. This helps
reduce the dimensionality of the explanatory variable space when categorical
explanatory variables need to be dummified in the modeling process.

  * •

For each of the remaining categorical explanatory variables, categories that
have too few observations are combined. For example, for the explanatory
variable STATE, which has the state information of the victim company, all
states that are associated with fewer than 1,000 observations are combined
into one category, named “Other”. This procedure helps mitigate the possible
problem that those small categories become absent in training set after the
train-test split. Moreover, with fewer categories, the dimensionality of the
sample is more manageable after the categorical explanatory variables are
dummified.

  * •

Auxiliary information, such as the IDs of companies and incidents, provides no
predictive power and thus is removed.

As a result, there are 8 explanatory variables remaining. Table
LABEL:tab:summary_stat shows the summary statistics of all these explanatory
variables; their descriptions are provided in Appendix A.

To provide an overview, the processed and cleaned sample contains 84,938
observations of cyber incidents, with the cyber incident type as the multi-
class categorical response variable, and 8 explanatory variables, of which 3
are numerical, and 5 are categorical. For categorical explanatory variables,
since many of them, such as state and industry, have more than two levels and
are non-ordinal, they are numerically coded by dummification. In addition, the
three numerical explanatory variables are highly right-skewed, as suggested by
their summary statistics, which are expected to have a negative impact on the
performance of some linear models that will be tested. Therefore, each of the
numerical explanatory variables is log-transformed. For testing and comparing
the performance of different models, 30% of all observations are selected at
random and used as the holdout sample, while the remaining 70% are for
training classifiers.

#### Multi-class classification

Predicting which type an incident falls into once it occurs is essentially a
multi-class classification problem. In this study, we explored several
commonly used classification models, including decision trees, random forests,
gradient-boosted trees, linear discriminant analysis, multinomial logistic
regression, and multi-layer perceptron. In the end, we build a stacking
classifier on top of these models for any possible performance gains. Because
all these methods are well established and documented in a broad volume of
literature, we shall only provide brief and mostly qualitative overviews of
them in Appendix B and focus instead on their performance and prediction
results.

#### Model comparisons

Table 1 suggests that the four incident categories are unbalanced, and a large
part of them are within the types of data breach and privacy violation types.
In this circumstance, several metrics, as summarized in [12], could be used to
evaluate model performance. In this paper, we shall report the balanced
accuracy of each model as a proxy of its performance. The balanced accuracy
score is defined as follows:

$\displaystyle\text{Balanced
Accuracy}=\frac{1}{2K}\Bigg{(}\sum_{k=1}^{K}\Bigg{(}$
$\displaystyle\;\frac{\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=\hat{O}_{i}=k\\}}}{\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=\hat{O}_{i}=k\\}}+\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=k\neq\hat{O}_{i}\\}}}$
(8) $\displaystyle+\frac{\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=\hat{O}_{i}\neq
k\\}}}{\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=\hat{O}_{i}\neq
k\\}}+\sum_{i=1}^{N}\mathbb{I}_{\\{\hat{O}_{i}=k\neq
O_{i}\\}}}\Bigg{)}\Bigg{)},$

where

  * •

$N$ is the size of the test set;

  * •

$O_{i}$ and $\hat{O}_{i}$ are the observed and the predicted occurring
incident types of record $i$, respectively, for $i=1,2,\dots,N$.

The one-versus-all metric evaluates the model performance by comparing the
predictions in one class against those in the combination of all other
classes. Specifically, for each $k\in\mathcal{K}$,
$\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=\hat{O}_{i}=k\\}}$ is the number of true
positives with respect to incident type $k$,
$\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=k\neq\hat{O}_{i}\\}}$ is the false
negative count, and $\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=\hat{O}_{i}\neq k\\}}$
and $\sum_{i=1}^{N}\mathbb{I}_{\\{\hat{O}_{i}=k\neq O_{i}\\}}$ correspond to
true negatives and false positives, respectively. Therefore,
$\frac{1}{2}\Bigg{(}\frac{\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=\hat{O}_{i}=k\\}}}{\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=\hat{O}_{i}=k\\}}+\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=k\neq\hat{O}_{i}\\}}}+\frac{\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=\hat{O}_{i}\neq
k\\}}}{\sum_{i=1}^{N}\mathbb{I}_{\\{O_{i}=\hat{O}_{i}\neq
k\\}}+\sum_{i=1}^{N}\mathbb{I}_{\\{\hat{O}_{i}=k\neq O_{i}\\}}}\Bigg{)}$
represents the $k$-th class-specific balanced accuracy, and Equation (8)
represents the average of these balanced accuracy values across all classes.

After all the classifiers are tuned and trained, their performance in terms of
the balanced accuracy of the predictions in the test set (holdout sample) is
presented in Table 2. The table shows the by-class balanced accuracy of each
classifier as well as the average across all classes, and the best score is
highlighted in each class or among all averages.

Classifier | By-Class Balanced Accuracy | Average
---|---|---
PV | DB | FE | ITE
Decision tree | 0.8373 | 0.7827 | 0.6369 | 0.6357 | 0.7231
Random forest | 0.8506 | 0.7976 | 0.7151 | 0.7193 | 0.7707
Gradient boosted trees | 0.8534 | 0.7993 | 0.7270 | 0.6991 | 0.7697
Linear discriminant analysis | 0.7005 | 0.6777 | 0.5244 | 0.5746 | 0.6193
Multinomial logistic regression | 0.7054 | 0.6817 | 0.4744 | 0.4864 | 0.5870
Multi-layer perceptron | 0.7515 | 0.7138 | 0.8571 | 0.4864 | 0.7022
Stack (Trees) | 0.8592 | 0.8021 | 0.7823 | 0.7799 | 0.8059
Stack (Linear classifiers excluded) | 0.8622 | 0.8030 | 0.7865 | 0.7698 | 0.8054
Stack (All) | 0.8616 | 0.8032 | 0.7760 | 0.7789 | 0.8049
Table 2: Comparison of classification models.

To summarize the performance comparison among different models, it is easy to
observe that stacking models overall outperform their individual base
classifiers. The best by-class accuracy scores and the best average accuracy
score are all achieved by stacking models. Second, the distinction between the
three stacking models is small. This suggests that the inclusion of
classifiers other than tree-based models does not substantially improve the
model performance. In addition, the by-class accuracy scores and the average
balanced accuracy score attained by stacking models are generally close to
$80\%$, which is a reasonably large proportion. Based on these three
observations, the stacking classifier built on top of the three tree-based
classifiers is adopted to predict incident-type occurrence probabilities.

The numerical results, in terms of the predicted probabilities of the four
incident types given different sets of company and incident characteristics,
are presented in Table 9 of Section 5.

###  3.3 Incident-specific loss severity

After predicting the occurrence probability of individual incident type, the
next step is to model the loss severity conditioning on a given incident type.
Therefore, another useful attribute of the dataset is the loss information on
each incident. Excluding observations that belong to the “other” category,
among all 84,373 incidents that belong to one of the 4 cyber incident
categories, there are 3,978 observations with known losses. The summary
statistics of those losses are shown in Table 3.

Type | Count | Min. | 1st Qu. | Median | Mean | 3rd Qu. | Max.
---|---|---|---|---|---|---|---
PV | 2,016 | 30 | 4,838 | 28,572 | 4,606,605 | 501,556 | 1,000,000,000
DB | 610 | 1 | 18,195 | 202,500 | 11,835,168 | 1,677,500 | 4,000,000,000
FE | 1,137 | 180 | 10,025 | 137,817 | 11,884,945 | 1,700,000 | 1,750,000,000
ITE | 215 | 200 | 20,000 | 194,850 | 4,994,540 | 1,050,000 | 200,000,000
All | 3,978 | 1 | 6,511 | 60,000 | 7,112,494 | 945,287 | 4,000,000,000
Table 3: Summary statistics of losses by incident type.

Prior to fitting any reasonable distributions to each conditional loss, it is
essential to realize that the conditional loss distributions could be
significantly different among cyber incident types. Table 4 shows the pairwise
two-sample Kolmogorov-Smirnov (KS) tests on empirical loss distributions among
various incident types. The numbers below the diagonal are the $p$-values of
these KS tests. Based on the significance level $0.05$, a $p$-value lower than
this level suggests that the null hypothesis, that the two empirical samples
being compared are from the same underlying distribution, should be rejected.
Due to symmetry, the statistical decisions are noted above the diagonal. The
result shows that the difference among the conditional losses of the four
incident types is statistically significant, and hence there is the necessity
of individually modeling each incident type’s loss severity, justifying the
incident-specific insurance coverage.

|  PV |  DB |  FE |  ITE
---|---|---|---|---
PV |  - |  reject |  reject |  reject
DB |  0.00E+00 |  - |  reject |  reject
FE |  0.00E+00 |  0.00E+00 |  - |  reject
ITE |  2.83E-05 |  1.70E-06 |  2.84E-05 |  -
Table 4: Pairwise comparison between empirical distributions of different
incident types using pairwise two-sample Kolmogorov-Smirnov test.

In addition, to show that loss severity distributions are mainly conditional
on incident types but not on additional explanatory variables, we performed an
ANOVA test with two linear regression models. One is a full model with total
loss as the response variable, and with incident type, as well as all
explanatory variables summarized in Table LABEL:tab:summary_stat. The other
model is a simple linear regression, and it has the total loss as the response
variable, but has incident type as its only explanatory variable. The ANOVA
test result, presented in Table 5, shows that the $p$-value of the test,
$0.1447$, is higher than any common choices of significance level, e.g.,
$0.05$, $0.01$, etc. Therefore, we can conclude that the additional
explanatory variables in Table LABEL:tab:summary_stat for the full model
contribute little to the explanation of the variance in losses. This lack of
explanatory power could be an artifact of the limited data size, which makes
it infeasible to build company-specific severity models, but this situation
can potentially be improved by a growing number of cyber incident records.

| Full | Simple
---|---|---
Residual Degrees of Freedom | $2783$ | $2829$
Residual Sum of Squares | $2.33\times 10^{19}$ | $2.38\times 10^{19}$
Difference in Degrees of Freedom |  | $-46$
Difference in Sum of Squares |  | $-4.71\times 10^{17}$
F Statistic |  | $1.2238$
$p$-value |  | $0.1447$
Table 5: ANOVA test between the full model and the model with only incident
type as its explanatory variable.

Figure 2 shows that the majority of the cases realize moderately small losses,
and a small number of cases realize substantial losses, indicating that the
conditional loss distributions of all incident types are likely to be right-
skewed. Therefore, we explore well-known distributions that could capture such
a skewness, namely log-normal, exponential, gamma, and Weibull, for each cyber
incident type.

(a) All losses.

(b) Losses greater than or equal to $10 million.

Figure 2: Histogram of losses of different incident types.

The parameters of each distribution are fitted using the maximum likelihood
estimation method. All fitted distributions for the same cyber incident type
are then compared based on their Akaike Information Criterion (AIC), which is
defined as

$\text{AIC}=2\hat{N}-2\ln(\hat{L}),$

where $\hat{N}$ is the number of parameters to be estimated, and $\hat{L}$ is
the maximum likelihood of the fitted model; a model with a smaller AIC value
is better. Table 6 summarizes the best-fitted loss distributions for all cyber
incident types. A comparison of the AICs of different fitted distributions is
presented in Table 12 of Appendix D.1. The distribution fitting results
suggest that it is most appropriate to treat losses of different incident
types as log-normal random variables with different log-mean and log-standard-
deviation parameters, as given in Table 6.

Type | Distribution | $\mu$ | $\sigma$
---|---|---|---
PV | log-normal | -2.5996 | 3.2798
DB | log-normal | -0.7916 | 3.1122
FE | log-normal | -3.4100 | 2.8577
ITE | log-normal | -1.9557 | 3.3629
Table 6: Best fitted loss distribution of each incident type.

###  3.4 Section summary

Thus far, all the necessary model inputs for solving Problem (5) are obtained
based on a dataset of historical cyber incidents. A brief recapitulation of
the key findings in this section is as follows:

  * •

In this study, the incident-specific cyber insurance policy would provide
coverage for four types of incidents, including privacy violation, data
breach, cyber fraud and extortion, and IT errors, i.e., $K=4$.

  * •

The occurrence probabilities of individual incident types are specific to each
company and incident, and the probability vector
$\bm{p}=(p_{1},p_{2},\dots,p_{K})$ is obtained by feeding the selected company
and incident features into a trained stacking classifier. This classifier
achieves a balanced accuracy score of around 80%, as defined in Equation (8),
for each incident type, which is superior to all other predictive models
surveyed (see Table 2).

  * •

Incident-specific losses are statistically significant to be different in
terms of severity. Additional explanatory variables, apart from the incident
type, do not seem crucial for modeling the loss severity. Log-normal
distributions provide the best fit for all conditional loss severity with
estimated parameters presented in Table 6.

This section demonstrates the feasibility of obtaining the model inputs
necessary to solve Problem (5) from real-world data. With approximately 80%
balanced accuracy scores associated with predictions of incident types and
well-fitted condition loss severity distributions, they are expected to be
practical and reliable to proceed to the next step.

##  4 Solution by Cross-Entropy Method

As discussed before, even with the necessary model inputs to solve Problem
(5), despite the minimization problem being finite-dimensional, it is also
combinatorial and, therefore, cannot be solved analytically to derive the
Pareto optimal $\bm{\theta}$ and $\bm{d}$, which are denoted as
$\bm{\theta}^{*}$ and $\bm{d}^{*}$. In this paper, we resort to the cross-
entropy method (CEM) proposed in [24, 25]. The following subsection provides a
brief review of the CEM adopting our own minimization problem in (5).

###  4.1 Brief review of cross-entropy method

To simplify the notations, let
$\mathbf{z}=\left(\bm{d},\bm{\theta}\right)\in\mathcal{Z}=\mathbb{R}_{+}^{4}\times\left\\{0,1\right\\}^{4}$,
and simply write the objective function in Problem (5) as
$F\left(\mathbf{z}\right)$ (keeping in mind that it depends on the
(conditional) distributions of $X$ and $O$ which are model inputs from Section
3). Thus, Problem (5) is
$\min_{\mathbf{z}\in\mathcal{Z}}F\left(\mathbf{z}\right)$. It is well-known
that an $\mathbf{z}^{*}\in\mathcal{Z}$ is optimal, if and only if, for any
$\varepsilon>0$, there exists an $\mathbf{z}\in\mathcal{Z}$ such that
$F\left(\mathbf{z}\right)<F\left(\mathbf{z}^{*}\right)+\varepsilon$; this is
also equivalent with requiring that for any $\varepsilon>0$,
$\mathbf{z}^{*}\in\mathcal{Z}_{\varepsilon}\left(\mathbf{z}^{*}\right)$, where
$\mathcal{Z}_{\varepsilon}\left(\mathbf{z}^{*}\right)=\left\\{\mathbf{z}\in\mathcal{Z}:F\left(\mathbf{z}\right)<F\left(\mathbf{z}^{*}\right)+\varepsilon\right\\}$.

The CEM is a numerical method to solve Problem (5) by adopting a probabilistic
approach to the characterization for optimality above. To this end, let
$\mathbf{Z}$ be a random vector with range $\mathcal{Z}$. The optimality can
then be further characterized by the fact that considering all possible
distributions of $\mathbf{Z}$, for any $\varepsilon>0$,
$\mathbb{P}\left(\mathbf{Z}\in\mathcal{Z}_{\varepsilon}\left(\mathbf{z}^{*}\right)\right)=\mathbb{P}\left(F\left(\mathbf{Z}\right)<F\left(\mathbf{z}^{*}\right)+\varepsilon\right)$
is maximized by the Dirac delta distribution $\delta_{\mathbf{z}^{*}}$ with
the optimal objective value being $1$ uniformly for all $\varepsilon>0$.
However, to maximize this associated probabilistic objective,
$\mathbf{z}^{*}$, the unknown to be solved for, has to be known ex-ante.
Hence, together with considering only parametric distributions of
$\mathbf{Z}$, the CEM solves the approximation counterpart of $\mathbf{z}^{*}$
as follows: for any $\tau$ being larger but close enough to
$F\left(\mathbf{z}^{*}\right)$, solve an
$\mathbf{u}^{*}\left(\tau\right)\in\mathcal{U}$ which maximizes
$\mathbb{P}\left(F\left(\mathbf{Z}\right)\leq\tau;\mathbf{u}\right)$, where
$\mathbf{Z}$ follows a distribution with a parametric density function
$f\left(\mathbf{z};\mathbf{u}\right)$, for $\mathbf{z}\in\mathcal{Z}$ and
$\mathbf{u}\in\mathcal{U}$, and $\mathcal{U}$ is the parameter set for the
distributions of $\mathbf{Z}$ under consideration. But still, without knowing
$\mathbf{z}^{*}$, the level $\tau$ in this approximation counterpart could not
be chosen appropriately.

The CEM proposes a two-phase approach to construct a sequence
$\left\\{\left(\hat{\tau}_{t},\hat{\mathbf{u}}_{t}\right)\right\\}_{t=1}^{\infty}$
such that it will converge to an $\left(\hat{\tau},\hat{\mathbf{u}}\right)$,
where $\left(\hat{\tau},\hat{\mathbf{u}}\right)$ is close to
$\left(F\left(\mathbf{z}^{*}\right),\mathbf{u}^{*}\right)$, and thus serves as
for the unknown $\left(F\left(\mathbf{z}^{*}\right),\mathbf{u}^{*}\right)$,
and where $\mathbf{u}^{*}\in\mathcal{U}$ such that
$f\left(\cdot;\mathbf{u}^{*}\right)$ is the corresponding Dirac density of the
Dirac delta distribution $\delta_{\mathbf{z}^{*}}$. The main steps of the two-
phase construction are as follows. Initialize an
$\hat{\mathbf{u}}_{0}\in\mathcal{U}$, and let $t$ iterate through $1,2,\dots$
with $\left(\hat{\tau}_{t},\hat{\mathbf{u}}_{t}\right)$ being updated in each
iteration. The $\hat{\tau}_{t}$, for $t=1,2,\dots$, is defined as the
$\varrho$-th quantile of the sample drawn from the distribution
$f\left(\cdot;\hat{\mathbf{u}}_{t-1}\right)$, where
$\varrho\in\left(0,1\right)$ is close to $0$ but is moderately small and
uniform for all $t=1,2,\dots$; in turn, the $\hat{\mathbf{u}}_{t}$, for
$t=1,2,\dots$, is defined as the maximizer of the following problem:

$\max_{\mathbf{u}_{t}\in\mathcal{U}}\hat{\mathbb{E}}\left[\mathbb{I}_{\left\\{F\left(\mathbf{Z}\right)\leq\hat{\tau}_{t}\right\\}}\ln
f\left(\mathbf{Z};\mathbf{u}_{t}\right);\hat{\mathbf{u}}_{t-1}\right],$ (9)

where $\hat{\mathbb{E}}\left[\cdot;\hat{\mathbf{u}}_{t-1}\right]$ represents
the estimated mean of the same sample drawn from the distribution
$f\left(\cdot;\hat{\mathbf{u}}_{t-1}\right)$. This two-phase update repeats
itself until the estimated variance of the next sample drawn from the updated
distribution $f\left(\cdot;\hat{\mathbf{u}}_{t}\right)$ is small enough and
the $\hat{\tau}_{t}$ does not change from those in the last few update steps
(see [5]). Therefore, the main idea of the construction is, first estimating
an unknown level $\hat{\tau}_{t}$ being larger but close enough to
$F\left(\mathbf{z}^{*}\right)$ as it is the $\varrho$-th sample quantile with
$\varrho$ being close to $0$ using the parameter $\hat{\mathbf{u}}_{t-1}$, and
second, estimating a parameter $\hat{\mathbf{u}}_{t}$ such that the event
$\left\\{F\left(\mathbf{Z}\right)\leq\hat{\tau}_{t}\right\\}$ with
$\mathbf{Z}\sim f\left(\cdot;\hat{\mathbf{u}}_{t}\right)$ is much more likely
to happen, compared to the case that $\mathbf{Z}\sim
f\left(\cdot;\hat{\mathbf{u}}_{t-1}\right)$ when its likelihood is given by a
moderately small $\varrho$, so to sequentially maximize
$\mathbb{P}\left(F\left(\mathbf{Z}\right)\leq\hat{\tau}_{t};\mathbf{u}_{t}\right)$.
In the second phase, the estimated parameter $\hat{\mathbf{u}}_{t}$ is to
minimize the cross-entropy between the unknown optimal density, for importance
sampling, and a reference density in the same parametric family; this is
equivalent to maximizing (9).

Note that in solving (9), only the observations in the sample drawn from the
distribution $f\left(\cdot;\hat{\mathbf{u}}_{t-1}\right)$ satisfying the event
condition, $F\left(\mathbf{Z}\right)\leq\hat{\tau}_{t}$, are useful to
calculate its objective; this motivates the introduction of an elite sample
consisting of observations in the original sample satisfying the condition.
Algorithm 1 provides the pseudo-code of the CEM.

Input: $\varrho\in(0,1)$ //Elite sample proportion

1

$f\left(\cdot;\cdot\right)$ //Parametric distribution of $\mathbf{Z}$

2

$\hat{\mathbf{u}}_{0}$ //Initial parameter for $f\left(\cdot;\cdot\right)$

3

$\nu$ //Variance stopping criterion

4

$l$ //Lag between iterations in which approximated optimums are compared

$N$ //Sample size

$F$ //Objective function

5

Output: Sample from $f\left(\cdot;\hat{\mathbf{u}}\right)$ //Approximation of
$\mathbf{z}^{*}$

6

7 $t\leftarrow 0$

$N_{E}\leftarrow\lceil\varrho N\rceil$ //Elite sample size

8

9Draw random sample
$\left\\{\mathbf{z}_{1}^{\left(t\right)},\mathbf{z}_{2}^{\left(t\right)},\dots,\mathbf{z}_{N}^{\left(t\right)}\right\\}$
from distribution with density $f(\cdot;\hat{\mathbf{u}}_{t})$

10Calculate element-wise sample variances $V^{\left(t\right)}_{k}$, for
$k=1,2,\dots,2K$

11$\Delta\leftarrow\infty$

12while
_$\max\left(V^{\left(t\right)}_{1},V^{\left(t\right)}_{2},\dots,V^{\left(t\right)}_{2K}\right)
>\nu$ or $\Delta\neq 0$ _ do

13 $y^{\left(t\right)}_{i}\leftarrow F(\mathbf{z}_{i}^{\left(t\right)})$ for
$i=1,2,\dots,N$

14 Select elite sample
$\mathcal{E}^{\left(t\right)}=\left\\{\mathbf{z}_{i}^{\left(t\right)}:F\left(\mathbf{z}_{i}^{\left(t\right)}\right)=y_{(i)}^{\left(t\right)},\;i=1,2,\dots,N_{E}\right\\}$

15 $\hat{\tau}_{t}\leftarrow y^{\left(t\right)}_{(N_{E})}$

16 if  _$t\geq l$_ then

17
$\Delta\leftarrow\max(|\hat{\tau}_{t}-\hat{\tau}_{t-1}|,|\hat{\tau}_{t}-\hat{\tau}_{t-2}|,\dots,|\hat{\tau}_{t}-\hat{\tau}_{t-l}|)$

18 end if

19

20 Obtain an estimate $\hat{\mathbf{u}}_{t+1}$ using the elite sample
$\mathcal{E}^{\left(t\right)}$

21 $t\leftarrow t+1$

22 Draw random sample
$\left\\{\mathbf{z}_{1}^{\left(t\right)},\mathbf{z}_{2}^{\left(t\right)},\dots,\mathbf{z}_{N}^{\left(t\right)}\right\\}$
from distribution with density $f(\cdot;\hat{\mathbf{u}}_{t})$

23 end while

24

25$T\leftarrow t$

return $\mathbf{z}_{1}^{\left(T\right)}$ //Arbitrarily choosing the first
observation because observations are similar due to small variance

26

Algorithm 1 Cross-Entropy Method

###  4.2 Implementation considerations and results

The actual implementation of this optimization problem consists of two parts,
including computing the objective function and implementing the CEM.

The objective function is the sum of the VaRs of both parties. As suggested by
Equations (6) and (7), the distribution of either the buyer’s or the seller’s
loss random variable is a mixture of distributions of incident-specific
losses, and the mixture distribution has no closed-form quantile function that
returns the Value-at-Risk of concern easily. Therefore, relying on the mixture
distribution function, we adopted a lookup approach to numerically find a
point that satisfies Equation (3) on a fine grid. In this study, risk
tolerance levels of the seller and the buyer are set at $\alpha=0.95$ and
$\beta=0.9$, respectively.

Although the CEM has a guaranteed asymptotic convergence (see, for example,
[26, 18]), the number of iterations needed to meet the stopping criteria is
uncertain because of the randomness introduced during the sampling processes.
Therefore, a common approach to addressing this issue in algorithm
implementation is adding additional stopping rules to ensure that the running
time of optimization tasks is manageable. Typical rules may include enforcing
a maximum number of iterations and assuming convergence if the value of the
objective function has not experienced a large enough improvement for a
predefined number of iterations. These rules, along with the randomness in
sampling, cause the algorithm in our implementation not always converge timely
and not always stop at the same value. As suggested by [5], the optimization
process shall be repeated several times for quality investigation and
assurance. Hence, for each set of incident probabilities predicted based on
certain company and incident characteristics, we ran 50 trials of the CEM to
solve the optimization with different random seeds, while holding other
specifications constant (see Appendix C). Results from all trials are stored
for further analysis.

Table 8 shows the results of five out of those 50 trials. Trial 5 converged to
an “optimum” larger than that achieved by Trial 1, suggesting that the
algorithm could possibly stop at a non-optimal point. This highlights the
necessity of multiple trials. Another observation is that, although most of
the trials reach a converging state after a few numbers of iterations (see,
e.g., Trials 1, 4, and 5), there exist initial random states that cause the
program to fail to meet the convergence criteria when the preset maximum
number of iterations is reached, e.g., Trials 2 and 3. Unless reducing the
maximum number of iterations at the risk of more non-convergent trials, the
great gap between numbers of iterations needed by different trials brings the
challenge that even in a multiprocessing environment, congestion occurs when
computational resources are occupied by those long-running trials, making
solving the optimization problem with repeating trials time-consuming.

Table 7 shows the amount of time spent on finding optimal insurance designs
for $5$ individual organizations using the CEM. The computation is done on 4
IBM POWER9 CPU cores. For each company, it takes around 2 minutes to complete
each trial. Given that we ran $50$ trials for each company, the running time
generally could easily exceed an hour.

| Organization
---|---
| 1 | 2 | 3 | 4 | 5
Count | 50 | 50 | 50 | 50 | 50
Mean (s) | 211.13 | 108.06 | 104.13 | 123.65 | 125.21
Standard Deviation (s) | 317.31 | 160.15 | 162.98 | 181.67 | 185.44
Minimum (s) | 24.97 | 11.68 | 9.62 | 15.41 | 14.84
25% Quantile (s) | 43.58 | 20.27 | 22.62 | 27.85 | 25.94
50% Quantile (s) | 68.18 | 34.78 | 36.63 | 45.37 | 40.52
75% Quantile (s) | 155.92 | 86.31 | 83.39 | 99.20 | 94.85
Maximum (s) | 975.61 | 538.80 | 539.01 | 612.82 | 542.56
Sum (s) | 10556.59 | 5402.97 | 5206.36 | 6182.30 | 6260.34
Table 7: Descriptive statistics on the number of seconds (s) spent on finding
the optimal insurance design using the CEM

Other than the long running time problem, the solution to the minimization
problem in Equation (5) is not unique, as presented in [4] for the case of two
incident types. The value of each element of $\bm{d}^{*}$ could lie in some
intervals, composed of points that result in the same minimum. Table 8 shows
that although both Trial 1 and 4 reached the converging state and attained the
same minimum, their solutions, $\bm{d}^{*}$ in particular, differ.
Nevertheless, the results provide common ground to insurers and policyholders
on the design of an incident-specific cyber insurance contract. As shown in
Table 8, based on the company’s own characteristics and the risk preferences
of the company and the insurer, both parties should agree that in
comprehensive cyber insurance covering PV, DB, FE, and ITE, a policy limit
should be implemented for PV, whereas deductibles should be applied on other
incident-specific coverages, since $\bm{\theta}^{*}=\left(0,1,1,1\right)$.
While holding the risk of both parties unchanged, deductible and limit amounts
can be negotiated for other possible considerations, one of which is discussed
in the next section.

Despite the computational challenges and non-unique solutions, the benefit of
such an incident-specific insurance contract is readily seen from the
comparison between the risks taken by both parties with and without insurance.
In the folloiwng, refer to Trials $1$ to $4$. If the company chooses not to
buy an incident-specific policy specified in Table 8, it will take the risk of
$13.6984 million. If it chooses to get covered, the appropriate amount of
premium $\pi$ (in millions) paid by the buyer should satisfy
$\$2.2977\leq\pi\leq\$4.6595$ according to Inequality (2), and that will yield
a risk reduction of $4.6595 million, which is always greater than the paid
premium. From the insurer’s perspective, the insurance contract will bring a
risk of $2.2977 million, but that will be fully compensated by the collected
premium. As a result, both parties are mutually benefited on their own
objectives; their aggregate risk is reduced by $2.3618 million, from $13.6984
million to $11.3366 million, which could be effective if the premium is agreed
to be lower than the reduced aggregate risk in its range. Again, note that the
premium discussed in this paper is derived economically which could be
actuarially priced in the insurance practice.

Trial | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$p_{1}^{*}$ (PV) | 0.3383
$p_{2}^{*}$ (DB) | 0.5717
$p_{3}^{*}$ (FE) | 0.0700
$p_{4}^{*}$ (ITE) | 0.0200
Random seed | 0 | 8 | 13 | 6 | 2
No. of iterations | 11 | 401 | 401 | 41 | 54
| Convergence
---
status
| Variance
---
converged
| Not
---
converge
| Not
---
converge
| Variance
---
converge
| Variance
---
converged
Optimum (millions) | 11.3366 | 11.3381
$\theta_{1}^{*}$ (PV) | 0 | 0
$\theta_{2}^{*}$ (DB) | 1 | 1
$\theta_{3}^{*}$ (FE) | 1 | 1
$\theta_{4}^{*}$ (ITE) | 1 | 1
$d_{1}^{*}$ (millions, PV) | 0.0531 | 0.0474 | 0.1074 | 0.0511 | 0.0334
$d_{2}^{*}$ (millions, DB) | 0.1011 | 0.1153 | 0.0982 | 0.1276 | 0.0941
$d_{3}^{*}$ (millions, FE) | 0.1167 | 0.0210 | 0.0107 | 0.0979 | 0.0520
$d_{4}^{*}$ (millions, ITE) | 0.1151 | 0.0195 | 0.0435 | 0.0314 | 0.0585
| Buyer’s risk without
---
insurance (millions)
13.6984 | 13.6984
| Buyer’s risk with
---
insurance (millions)
9.0389 | 9.0389
| Buyer’s risk
---
reduction (millions)
4.6595 | 4.6595
| Seller’s risk without
---
insurance (millions)
0.0000 | 0.0000
| Seller’s risk with
---
insurance (millions)
2.2977 | 2.2992
| Seller’s risk
---
increase (millions)
2.2977 | 2.2992
| Aggregate risk without
---
insurance (millions)
13.6984 | 13.6984
| Aggregate risk with
---
insurance (millions)
11.3366 | 11.3381
| Aggregate risk
---
reduction (millions)
2.3618 | 2.3603
| Premium range
---
(millions)
[2.2977, 4.6595] | [2.2992, 4.6595]

Table 8: Results of five trials of CEM with the same set of predicted incident
probabilities.

##  5 Function Approximation

In Section 4.2, two computational challenges are discussed in the
implementation of the optimization algorithm, including the numerical
evaluation of quantiles of loss random variables and random states that lead
to non-convergent trials. Both challenges induce long running time and could
make this policy design less practical in the production environment. To speed
up the computation process, we make use of function approximation to find a
target function, which establishes a more direct mapping between incident
properties, including their occurrence probabilities and severities, and the
parameters of the eventual indemnity functions in the incident-specific
policy. Note that in this study, the distribution parameters of incident
severities are fixed, and only probabilities vary depending on company
characteristics; therefore, more precisely speaking, the mapping is only
between incident probabilities and indemnity function parameters.

In this section, we distinguish between true solutions, which are obtained by
using CEM to approximately solve Problem (5), and fitted solutions, which are
obtained through function approximation. By using fitted solutions to
approximate true solutions, such that the fitted solutions result in the same
optimum in Problem (5) as the true solutions, the end goal is that once a
company’s characteristics are fed into the trained stacking classifier built
in Section 3.2, the corresponding probability vector can be used as inputs of
the approximate function to generate the parameters of the optimal incident-
specific cyber insurance, i.e., $\hat{\bm{\theta}}$ and $\hat{\bm{d}}$.

###  5.1 Sample selection

Provided that the true solutions to Problem (5) are non-unique, as shown in
Section 4, the mapping between $\bm{p}$ and
$\left(\bm{\theta}^{*},\bm{d}^{*}\right)$ is one-to-many. This could
potentially lead to a lack of fit and suboptimal fitted solutions. The
suboptimality is a result of the fact that the true solutions of individual
$d^{*}_{k}$, for $k=1,\dots,K$, are on disjoint intervals on the real line
(see [4]) but the function approximation process can easily disregard that
premise if, with respect to the same set of probabilities, the divergence
between the fitted solution and all true solutions, e.g., mean squared error
or mean absolute error, is to be minimized because the fitted value could
plausibly be outside the intervals where the true solutions reside in.

Therefore, instead of looking for a fitted solution with the smallest average
deviation from all true solutions, it is more guaranteed to find an optimal
fitted solution close to one of the true solutions. In this study, to identify
a unique true solution, we choose the one that minimizes the seller’s expected
loss among non-unique true solutions. That is, the solution selected for
function approximation is, therefore,
${\arg\min}_{(\bm{d},\bm{\theta})\in\mathcal{S}}\mathbb{E}L_{s}(\bm{d},\bm{\theta};X,O)$,
where $\mathcal{S}$ is the set of Pareto optimal
$\left(\bm{\theta}^{*},\bm{d}^{*}\right)$ computed by the CEM. This strategy
serves as an illustration of how a unique true solution can be chosen based on
additional preferences, while in practice, alternative preferences can be used
to compare among the non-unique solutions to Problem (5) and choose the best
one.

###  5.2 Function approximation procedure

Because $\bm{d}$ depends on $\bm{\theta}$, the function approximation follows
a two-step process. The relationship between probabilities $\bm{p}$ and the
types of indemnities $\bm{\theta}$ is first approximated, for which $\bm{p}$
are seen as features, and $\bm{\theta}$ are seen as labels. Because each
$\theta_{k}$, for $k=1,\dots,K$, takes the value of either $0$ or $1$ and is
not mutually exclusive, predicting $\bm{\theta}$ is a multi-label
classification problem. Then, conditioning on $\bm{\theta}$, the mapping
between probabilities $\bm{p}$ and the amounts of incident-specific
deductibles or limits $\bm{d}$ shall be approximated, and this is a multi-
output regression problem. Because there are $2^{K}$ distinct values of
$\bm{\theta}$, the number of regression models to build is $2^{K}$. The target
function to be approximated then consists of both the classification model and
the conditional regression models.

Let $G\colon(0,1)^{K}\rightarrow\\{0,1\\}^{K}$ be the classifier and
$H_{\bm{\theta}}\colon(0,1)^{K}\rightarrow\mathbb{R}_{+}^{K}$ be the
conditional regression models depending on $\bm{\theta}$. In addition, the
sample created in Section 5.1 is split into a training set
$\\{(\bm{p}_{s};\bm{\theta}^{*}_{s},\bm{d}^{*}_{s})\\}_{s=1}^{S_{1}}$ and a
test set
$\\{(\bm{p}_{S_{1}+s};\bm{\theta}^{*}_{S_{1}+s},\bm{d}^{*}_{S_{1}+s})\\}_{s=1}^{S_{2}}$,
where $S_{1}$ and $S_{2}$ are the sizes of the two sets, respectively. The
function approximation procedure can be summarized in Algorithm 2.

//Training

Input: $\\{(\bm{p}_{s};\bm{\theta}^{*}_{s},\bm{d}^{*}_{s})\\}_{s=1}^{S_{1}}$,
initial classifier $G^{0}$, initial regression models $H^{0}_{\bm{\theta}}$

1

Output: Trained classifier $G^{*}$, trained regression models
$H^{*}_{\bm{\theta}}$

2

3Training the multi-label classifier $G^{0}$ on the training set with $\bm{p}$
as features and $\bm{\theta}$ as labels

4for  _ $\bm{\theta}\in\\{0,1\\}^{K}$_ do

5 Training the multi-output regression model $H^{0}_{\bm{\theta}}$ on the
training set with $\bm{p}$ as features and $\bm{d}$ as labels

6 end for

7

//Testing

Input:
$\\{(\bm{p}_{S_{1}+s};\bm{\theta}_{S_{1}+s}^{*},\bm{d}_{S_{1}+s}^{*})\\}_{s=1}^{S_{2}}$,
trained classifier $G^{*}$, trained regression models $H^{*}_{\bm{\theta}}$

8

Output: Error rate, fitted solutions
$\\{(\hat{\bm{\theta}}_{S_{1}+s},\hat{\bm{d}}_{S_{1}+s})\\}_{s=1}^{S_{2}}$

9

10for  _$s=1,2,\dots,S_{2}$_ do

11 $\hat{\bm{\theta}}_{S_{1}+s}\leftarrow G^{*}(\bm{p}_{S_{1}+s})$

12 $\hat{\bm{d}}_{S_{1}+s}\leftarrow
H_{\hat{\bm{\theta}}_{S_{1}+s}}^{*}(\bm{p}_{S_{1}+s})$

13 $\epsilon_{S_{1}+s}\leftarrow
F(\bm{d}^{*}_{S_{1}+s},\bm{\theta}^{*}_{S_{1}+s};X,O)-F(\hat{\bm{d}}_{S_{1}+s},\hat{\bm{\theta}}_{S_{1}+s};X,O)$

14 end for

15

Error rate
$\leftarrow\left(\sum_{s=1}^{S_{2}}\mathbb{I}_{\\{\epsilon_{S_{1}+s}\neq
0\\}}\right)\big{/}{S_{2}}$

Algorithm 2 Function Approximation: Model training and testing

For the evaluation of the fitted solutions produced on the test set, error is
defined as the difference between the values of the objective function in
Problem (5) when true solutions and fitted solutions are provided,
respectively, i.e.,
${\epsilon_{s}=F(\bm{d}^{*}_{s},\bm{\theta}^{*}_{s};X,O)-F(\hat{\bm{d}}_{s},\hat{\bm{\theta}}_{s};X,O)}$.
The error rate is used as a performance measure of the target function, which
represents the ratio of predictions made on the test set resulting in the same
optimums as true solutions.

###  5.3 Numerical results

For each distinct set of probabilities $\bm{p}$, $50$ corresponding true
solutions are computed using the CEM, each of which has a different random
state, and then among them, a unique true solution, which minimizes the
seller’s expected loss, is selected for the sample in function approximation.

The computation of individual optimization tasks with specifications as shown
in Appendix C is parallelized on $96$ CPU cores. During each 24-hour wall-
clock time, the average number of unique true solutions that can be generated
is around $790$. This computationally intense process highlights the necessity
of the use of function approximation in practice. For both the classification
and regression tasks, we build tree-based models using the LightGBM framework
(see [16]), with a training set of size 2,201. Table 9 shows the fitted
solutions of five organizations, provided their organizational
characteristics.

Organization | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
NCASE | 55 | 652 | 28 | 3300 | 270
EMP | 4500 | 25105 | 119 | 16000 | 7908
REV (millions) | 461.47 | 40653.00 | 0.72 | 2662.68 | 1448.391
CTYPE | Private | Public | Private | Public | Other
STATE | Other | California | Arizona | Connecticut | Florida
IND | Services | Services | Retail Trade |  | Transportation &
---
Public Utilities
| Public
---
Administration
Incident Occurrence
AFTER2012 | Yes | Yes | Yes | Yes | Yes
MON | January | January | January | March | March
$p_{1}$ | 0.3383 | 0.4401 | 0.4700 | 0.4340 | 0.2300
$p_{2}$ | 0.5717 | 0.3340 | 0.3400 | 0.4360 | 0.4800
$p_{3}$ | 0.0700 | 0.1764 | 0.1600 | 0.0600 | 0.1900
$p_{4}$ | 0.0200 | 0.0495 | 0.0300 | 0.0700 | 0.1000
$\hat{\theta}_{1}$ | 0 | 1 | 1 | 0 | 0
$\hat{\theta}_{2}$ | 1 | 1 | 1 | 1 | 1
$\hat{\theta}_{3}$ | 1 | 0 | 0 | 1 | 1
$\hat{\theta}_{4}$ | 1 | 0 | 0 | 1 | 0
$\hat{d}_{1}$ (millions) | 0.1430 | 0.0968 | 0.0943 | 0.1318 | 0.1094
$\hat{d}_{2}$ (millions) | 0.1502 | 0.0901 | 0.0924 | 0.1636 | 0.1852
$\hat{d}_{3}$ (millions) | 0.0957 | 0.0904 | 0.0781 | 0.0918 | 0.1477
$\hat{d}_{4}$ (millions) | 0.0768 | 0.1490 | 0.1249 | 0.0974 | 0.1228
| CEM optimum
---
(millions)
11.3366 | 7.4182 | 7.7513 | 9.7930 | 9.2889
| Function
---
approximation
optimum (millions)
11.3366 | 7.4197 | 7.7513 | 9.7930 | 9.2889
| Error
---
(millions)
0.0000 | 0.0015 | 0.0000 | 0.0000 | 0.0000
| Buyer’s risk without
---
insurance (millions)
13.6984 | 8.5989 | 8.7049 | 11.4981 | 11.3691
| Buyer’s risk with
---
insurance (millions)
9.0389 | 6.9697 | 7.4837 | 5.9956 | 6.5337
| Buyer’s risk
---
reduction (millions)
4.6595 | 1.6292 | 1.2212 | 5.5025 | 4.8354
| Seller’s risk without
---
insurance (millions)
0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000
| Seller’s risk with
---
insurance (millions)
2.2977 | 0.4500 | 0.2676 | 3.7974 | 2.7552
| Seller’s risk
---
increase (millions)
2.2977 | 0.4500 | 0.2676 | 3.7974 | 2.7552
| Aggregate risk without
---
insurance (millions)
13.6984 | 8.5989 | 8.7049 | 11.4981 | 11.3691
| Aggregate risk with
---
insurance (millions)
11.3366 | 7.4197 | 7.7513 | 9.7930 | 9.2889
| Aggregate risk
---
reduction (millions)
2.3618 | 1.1792 | 0.9536 | 1.7051 | 2.0802
| Premium range
---
(millions)
[2.2977, 4.6595] | [0.4500, 1.6292] | [0.2676, 1.2212] | [3.7974, 5.5025] | [2.7552, 4.8354]

Table 9: Comparisons between exact solutions solved by the CEM and fitted
solutions by function approximation.

This table summarizes multiple aspects of the workflow proposed by this paper.
First, there are indeed different optimal incident-specific insurance designs
for organizations of distinct characteristics. The presented organizations
vary in litigation history, size, ownership, location, and industry. These
characteristics can be used as rating factors in the underwriting process, and
the insurer should be able to collect these pieces of information on the
insurance buyer conveniently. This result validates the practice of coverage
being designed in an incident-specific manner.

Second, similar to the results shown in Table 8 of Subsection 4.2, Table 9
demonstrates how both the insured and the insurer can be mutually benefited
from an incident-specific cyber policy. Indeed, Organization 1 is the same
organization used to derive the results in Table 8, and the numbers regarding
risks with and without insurance and the premium range are identical in both
tables. We shall relegate readers to Subsection 4.2 for a detailed account of
the policy being mutually beneficial to both parties. In addition, the same
conclusion can be drawn from the other four organizations presented in Table
9, suggesting that the proposed incident-specific policy can offer such a
benefit to organizations of different characteristics. However, Table 9 shows
that the effectiveness of aggregate cyber risk reduction does not always hold.
Indeed, while the premiums to be charged for Organizations 1, 2, and 3 could
be lower than their respective aggregate cyber risk reductions, the premiums
to be charged for Organizations 4 and 5 must be larger than their respective
aggregate risk reductions.

Lastly, in most cases, the target function can produce solutions that yield
the exact optimum as the true solutions generated by the CEM. A small
percentage of the fitted solutions are non-optimal; see, for example,
Organization 2. This problem can potentially be mitigated by employing a
larger training set to train the target function. Overall, an error rate of
$0.22$ on the test set of size 790 is attained. That is, 78% of the fitted
solutions result in the same optimal risks as those given by true solutions.
The slight compromise in accuracy leads to high computational efficiency. For
a test set of size $5$, the time needed for computing the fitted solutions is
approximately $0.068$ seconds. The running time here and the running time of
the CEM presented in Table 7 are measured on the same hardware. Compared to
the aforementioned CEM, which may take hours to solve the problem, the usage
of function approximation makes it much more feasible to generate policy
specifications on the fly in the production environment.

##  6 Conclusions and Future Directions

In this paper, we proposed a workflow for the design of incident-specific
cyber insurance. It consists of three key components, including the estimation
of model inputs using public and proprietary data in the underwriting process,
the Pareto optimal objective based on the insured’s and the insurer’s
preference orderings by their risk measures, and solvers for the optimization
problem. Using real cyber incident data, we showed how this workflow generates
incident-specific policies that lower the total risk taken by the insured and
the insurer. We also demonstrated that the proposed workflow can be time
efficient with the help of function approximation if there is a time
constraint on the delivery of results, which is common for insurers that
provide quotes and coverages on the fly.

Despite the time efficiency of using an approximated target function, one
limitation is that the optimum is not always attained by the fitted solutions
generated by the target function. Therefore, the benefits and risks of using
this approach should be weighed in practice. Another limitation of this study
is that, because of the scarcity of data on cyber incident losses, it is not
feasible to model company-specific loss distributions for individual incident
types, and therefore we only considered how incident type realizations are
related to company characteristics, regardless of the relationship between
those characteristics and the potential cyber losses. This issue could
possibly be resolved in future studies if a better quality of cyber loss
dataset is available. Lastly, this study assumes that indemnity functions are
specific to cyber incident types, or equivalently perils, only. In practice,
deductibles and limits could also be placed on more granular levels like
covered assets, such as hardware, data, or business income; see, for example,
[7]. Depending on the use case, the proposed workflow can easily be adopted
for policies that offer coverages at such more granular level.

## Acknowledgments

Wing Fung Chong, Daniël Linders, and Linfeng Zhang are supported by a 2018
Ignacio Hernando de Larramendi Research Grant from the Fundación MAPFRE. Wing
Fung Chong, Daniël Linders, Zhiyu Quan, and Linfeng Zhang are supported by a
Centers of Actuarial Excellence (CAE) Research Grant (2019–2021) from the
Society of Actuaries (SOA). Any opinions, findings, conclusions, or
recommendations expressed in this material are those of the authors and do not
necessarily reflect the views of the Fundación MAPFRE and the SOA.

## References

  * [1] AIG  “CyberEdge® Policy Documentation”, 2019  URL: https://www.aig.co.uk/content/dam/aig/emea/united-kingdom/documents/Financial-lines/Cyber/cyberedge-policy-documentation.pdf
  * [2] Zeinab Amin  “A Practical Road Map for Assessing Cyber Risk”  In _Journal of Risk Research_ 22.1 Informa UK Limited, 2017, pp. 32–43
  * [3] Alexandru V. Asimit and Tim J. Boonen  “Insurance with multiple insurers: A game-theoretic approach”  In _European Journal of Operational Research_ 267.2, 2018, pp. 778–790
  * [4] Alexandru V. Asimit, Tim J. Boonen, Yichun Chi and Wing Fung Chong  “Risk sharing with multiple indemnity environments”  In _European Journal of Operational Research_ 295.2, 2021, pp. 587–603
  * [5] Tim Benham, Qibin Duan, Dirk P. Kroese and Benoît Liquet  “CEoptim: Cross-Entropy R Package for Optimization”  In _Journal of Statistical Software_ 76.8 Foundation for Open Access Statistic, 2017
  * [6] James Cebula and Lisa Young  “A Taxonomy of Operational Cyber Security Risks”, 2010  URL: http://resources.sei.cmu.edu/library/asset-view.cfm?AssetID=9395
  * [7] Wing Fung Chong, Runhuan Feng, Hins Hu and Linfeng Zhang  “Cyber Risk Assessment for Capital Management”, 2022  arXiv:2205.08435 [q-fin.RM]
  * [8] Michel Dacorogna, Nehla Debbabi and Marie Kratz  “Building up cyber resilience by better grasping cyber risk via a new algorithm for modelling heavy-tailed data”  In _European Journal of Operational Research_ 311.2, 2023, pp. 708–729
  * [9] Martin Eling and Jan Wirfs  “What Are the Actual Costs of Cyber Risk Events?”  In _European Journal of Operational Research_ 272.3 Elsevier BV, 2019, pp. 1109–1119
  * [10] Equifax  “2017 Cybersecurity Incident & Important Consumer Information”, 2017  URL: https://www.equifaxsecurity2017.com/consumer-notice/
  * [11] Abhijeet Ghadge, Maximilian Weiß, Nigel D. Caldwell and Richard Wilding  “Managing cyber risk in supply chains: a review and research agenda”  In _Supply Chain Management: An International Journal_ 25.2 Emerald, 2019, pp. 223–240
  * [12] Margherita Grandini, Enrico Bagli and Giorgio Visani  “Metrics for Multi-Class Classification: an Overview”, 2020  arXiv:2008.05756 [stat.ML]
  * [13] Jay Greene  “Amazon’s cloud-computing outage on Wednesday was triggered by effort to boost system’s capacity”, 2020  URL: https://www.washingtonpost.com/technology/2020/11/28/amazon-outage-explained
  * [14] Jake Holland  “Facebook’s $650 Million Privacy Settlement Approved By Judge”, 2021  URL: https://news.bloomberglaw.com/privacy-and-data-security/650-million-facebook-privacy-settlement-approved-by-judge
  * [15] Tom Johansmeyer  “Cybersecurity Insurance Has a Big Problem”  In _Harvard Business Review_ , 2021  URL: https://hbr.org/2021/01/cybersecurity-insurance-has-a-big-problem
  * [16] Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye and Tie-Yan Liu  “Lightgbm: A highly efficient gradient boosting decision tree”  In _Advances in neural information processing systems_ 30, 2017, pp. 3146–3154
  * [17] Jay P. Kesan and Linfeng Zhang  “Analysis of Cyber Incident Categories Based on Losses”  In _ACM Transactions on Management Information Systems_ 11.4 Association for Computing Machinery (ACM), 2020, pp. 1–28
  * [18] L. Margolin  “On the Convergence of the Cross-Entropy Method”  In _Annals of Operations Research_ 134.1 Springer ScienceBusiness Media LLC, 2005, pp. 201–214
  * [19] Angelica Marotta, Fabio Martinelli, Stefano Nanni, Albina Orlando and Artsiom Yautsiukhin  “Cyber-insurance survey”  In _Computer Science Review_ 24 Elsevier BV, 2017, pp. 35–61
  * [20] Christopher J. Merz  “Using Correspondence Analysis to Combine Classifiers”  In _Machine Learning_ 36.1/2 Springer ScienceBusiness Media LLC, 1999, pp. 33–58
  * [21] NIST  “Framework for Improving Critical Infrastructure Cybersecurity”, 2018  National Institute of StandardsTechnology  URL: https://nvlpubs.nist.gov/nistpubs/CSWP/NIST.CSWP.04162018.pdf
  * [22] Office of Cybersecurity, Energy Security, and Emergency Response  “Colonial Pipeline Cyber Incident”, 2021  URL: https://www.energy.gov/ceser/colonial-pipeline-cyber-incident
  * [23] Sasha Romanosky, Lillian Ablon, Andreas Kuehn and Therese Jones  “Content Analysis of Cyber Insurance Policies: How Do Carriers Price Cyber Risk?”  In _Journal of Cybersecurity_ 5.1 Oxford University Press (OUP), 2019
  * [24] Reuven Rubinstein  “The Cross-entropy Method for Combinatorial and Continuous Optimization”  In _Methodology And Computing In Applied Probability_ 1.2 Springer ScienceBusiness Media LLC, 1999, pp. 127–190
  * [25] Reuven Y. Rubinstein  “Combinatorial Optimization, Cross-entropy, Ants and Rare Events”  In _Applied Optimization_ Springer US, 2001, pp. 303–363
  * [26] Reuven Y. Rubinstein  “Cross-entropy and rare events for maximal cut and partition problems”  In _ACM Transactions on Modeling and Computer Simulation_ 12.1 Association for Computing Machinery (ACM), 2002, pp. 27–53
  * [27] Daniel Woods, Ioannis Agrafiotis, Jason R.. Nurse and Sadie Creese  “Mapping the coverage of security controls in cyber insurance proposal forms”  In _Journal of Internet Services and Applications_ 8.1 Sociedade Brasileira de Computacao - SB, 2017

##  Appendix A Explanatory Variables and Their Descriptions

Explanatory Variable | Description
---|---
EMP | Number of employees
NCASE | Number of Federal Docket cases linked to company
REV | Total revenue (in millions of USD)
MON | Month in which the incident occurs
STATE | State in which the victim company is based
CTYPE | Company type
| * PRV (Private)
| * PUB (Public)
IND | Industy by SIC divisions
YEAR | Whether the incident occurred before 2012 or after 2012
Table 10: All explanatory variables and their descriptions.

##  Appendix B Classification Models

##### Decision tree

A decision tree (DT) splits the dataset based on the values of explanatory
variables and attempts to partition the dataset into homogeneous subsets. Each
split of a parent node is performed with the goal of reducing the Gini
impurity of data in child nodes. Decision tree is known to have a tendency to
overfit, which means that it can fit the training data well by growing a
complex tree with granular partitions, but its performance can suffer when
making predictions on new data. To mitigate this issue, three hyperparameters,
including the maximum depth of the tree, the minimum number of samples
required in each leaf, and the complexity parameter of the tree pruning
process, are tuned by grid search with 5-fold cross validations, to improve
the out-of-sample prediction performance of the built tree. The hyperparameter
tuning process of all other classification methods mentioned in this paper are
performed in the same way.

##### Random forest

Random Forest (RF) is an ensemble model that employs multiple classification
trees, of which each is built based on a bootstrapped sample of the training
set and a randomly selected subset of all explanatory variables. In the end,
the probabilistic predictions made by all trees are averaged to give a final
prediction. The randomness introduced by bootstrapping and randomly selected
explanatory variables help reduce the variance of the model, thus mitigating
the overfitting issue of a single classification tree, as discussed
previously. Hyperparameters including the number of trees to build and the
size of the subset of explanatory variables are tuned for performance
improvement in this study.

##### Gradient boosted trees

Gradient boosted trees is another type of ensemble model based on decision
trees. Instead of growing trees independently and then aggregating the
results, the boosted method iteratively improves the performance of the
previous tree-based classifier by updating it with a new tree that fits the
residual. The updating process is a linear combination of the old classifier
and the newly built tree, in which the coefficients minimize a certain loss
function of the true values and the predicted values. Hyperparameters,
including the learning rate, which governs how much each new tree contribute
to the final classifier, and the total number of trees to build, are tuned.

##### Linear discriminant analysis

Linear discriminant analysis (LDA) solves an optimization problem, for which a
transformation matrix is solved such that after the transformation, the ratio
between the total between-class variance of data points measured in an
Euclidean space and the total within-class variance is maximized. That is,
data points in the same class are made as close to each other as possible,
whereas data points belonging to different classes are separated to the
largest extent. The advantage of this method over the three tree-based models
introduced earlier is that there is an analytical solution to the
transformation matrix, and thus the training time is minimal. The drawback,
however, is that it imposes relatively strict assumptions on the data, such as
multivariate normality, independence, and non-collinearity, which can rarely
be fully satisfied by empirical data, and in that case, the performance can
suffer. Given that our sample contains many dummy variables, the assumptions
are largely violated. Nevertheless, this method is tested and compared with
other models.

##### Multinomial logistic regression

For multi-class problems, multinomial logistic regression (MLR) is a
generalization of its binary counterpart. Let $K$ denote the number of classes
and for each class $k=1,2,\dots,K$, the probability $p_{k}$ of an observation
with features vector $\bm{x}$ belonging to that class is modeled as follows,

$\ln p_{k}=\beta_{k}\cdot\bm{x}-\ln\sum_{k=1}^{K}e^{\beta_{k}\cdot\bm{x}},$

where $\beta_{k}$ is the coefficient learned from the training set to minimize
the cross-entropy loss function, and the term
$\ln\sum_{k=1}^{K}e^{\beta_{k}\cdot\bm{x}}$ is a normalizing factor that
ensures $\sum_{k=1}^{K}p_{k}=1$. Equivalently,

$p_{k}=\frac{e^{\beta_{k}\cdot\bm{x}}}{\sum_{k=1}^{K}e^{\beta_{k}\cdot\bm{x}}},\;\text{
for }k=1,2,\dots,K$

which is the softmax function of $\beta_{k}\cdot\bm{x}$.

##### Multi-layer perceptron

Multinomial logistic regression can be viewed as a single-layer perceptron
without hidden layers, of which the activation function in the output layer is
softmax, as introduced earlier. Thus, it would be interesting to compare
multinomial logistic regression to multi-layer perceptron (MLP) models with
additional hidden layers. Hidden layers enable the model to learn non-linear
relationships, and thus might provide a better fitting. For experiment
purposes, an MLP classifier with two hidden layers and ten neurons in each
hidden layer is trained.

##### Stacking classifier

Lastly, stacking classifiers are built based on the results output by the
different algorithms mentioned above. The motivation is that classifiers of
different kinds are likely to produce errors in different ways, and thus
combining their results with proper weights may result in reduced errors, as
suggested in [20]. To stack the base classifiers as introduced previously,
their probabilistic predictions are collected and used as new features for
building meta-classifiers on top of all individual models. The choice of the
meta-classifier in this study is random forest. To see how different sets of
base classifiers can affect the performance of the meta-classifier, three
stacking classifiers are compared. The first uses only tree-based base
classifiers, including decision trees, random forests, and gradient-boosted
trees. The second stacking classifier excludes the two linear methods, which
are linear discriminant analysis and multinomial logistic regression. The
third stacking classifier includes all models mentioned previously as base
classifiers.

##  Appendix C Optimization Parameters

Model Inputs |  Value
---|---
Seller’s Risk Level $\alpha$ |  0.95
Buyer’s Risk Level $\beta$ |  0.9
CEM Specifications |  Value
Initial Distribution of $\bm{d}$ |  Truncated multivariate normal with the mean and standard deviation of each marginal distribution to be 0 and 100,000 respectively
Initial Distribution of $\bm{\theta}$ |  Multivariate Bernoulli with the success probability of each marginal distribution to be 0.5
Variance threshold of $\bm{d}$ |  0.1
Variance threshold of $\bm{\theta}$ |  0.01
Sample size $N$ |  10
Elite sample proportion $\varrho$ |  0.2
No-improve iterations before stop $l$ |  10
Maximum Number of iterations |  401
Table 11: Seller’s and buyer’s risk levels and CEM specifications

##  Appendix D Severity Distribution Fitting Results

###  D.1 Comparison of AICs of fitted distributions

| Privacy Violation | Data Breach | Extortion/Fraud | IT Error
---|---|---|---|---
log-normal | 32.82 | 2154.44 | -2135.91 | 294.72
exponential | 9422.81 | 5014.21 | 4259.55 | 1714.78
gamma | 939.88 | 2484.55 | -1167.95 | 454.25
weibull | 269.59 | 2228.09 | -1896.31 | 352.69
Table 12: AICs of distributions fitted to incident-specific losses. For each
incident type, the log-normal distribution has the lowest value among all
fitted distributions.