Robust Actuarial Risk Analysis
Jose Blanchet∗

Henry Lam†

Qihe Tang‡

Zhongyi Yuan§

Abstract
This paper investigates techniques for the assessment of model error in the context of insurance risk analysis. The methodology is based on finding robust estimates for actuarial quantities
of interest, which are obtained by solving optimization problems over the unknown probabilistic
models, with constraints capturing potential nonparametric misspecification of the true model.
We demonstrate the solution techniques and the interpretations of these optimization problems,and illustrate several examples including calculating loss probabilities and conditional value-atrisk.

1

Introduction

This paper studies a methodology to quantify the impact of model assumptions that are uncertain
or possibly incorrect in actuarial risk analysis. The motivation of our study is that, in many
situations, coming up with a highly accurate model is challenging. This could be due to a lack ofdata, e.g., in describing a low-probability event, or modeling issues, e.g., in addressing a hidden
and potentially sophisticated dependence structure among risk factors. Often times, actuaries andstatisticians resort to models and calibration procedures that capture stylized features based on
experience or expert knowledge, but bearing the risk of deviating too much from reality and thus
leading to suboptimal decision-making. The methodology studied in this paper aims to quantify
the incurred errors from these approaches to an extent that we shall describe.
On a high level, the methodology we study in this paper has the following characteristics:∗

Department of Management Science and Engineering, Stanford University, Stanford, CA.
Department of Industrial Engineering and Operations Research, Columbia University, New York, NY.
‡
Department of Statistics and Actuarial Science, University of Iowa, Iowa City, IA.
§
Department of Risk Management, Pennsylvania State University, University Park, PA.
†121) Its starting point is a baseline model that, for any good reasons (e.g., a balance of fidelity
and tractability), is currently employed by the actuary.
2) The method employs optimization theory to find a bound for the underlying risk metric ofinterest, where the optimization is non-parametrically imposed over all probabilistic models thatare in the neighborhood of the baseline model.
3) Typically the worst-case model (i.e., the optimal solution obtained from the optimization
procedure) can be written in terms of the baseline model. So the proposed procedure can beunderstood as a correction to allow for the possibility of model misspecification.
The key idea of the methodology lies in a suitable formulation of the optimization problem in2) that leads to 3). This is an optimization cast over the space of models. This formulation has afeasible region described by a neighborhood of the baseline model, the latter believed to be close toa “true” or correct one - thanks to the expertise of the actuary. However, despite in the vicinity ofthe baseline, the location of the true model is unknown. The neighborhood that defines the feasible
region, roughly speaking, serves to provide a more accessible indication of where the true model is.
When this neighborhood happens to contain the true model, the optimization problem, namely amaximization or a minimization of the risk metric, will output an upper or a lower bound for thecorresponding true value.
Thus, instead of outputting a risk metric value obtained from a single baseline model, thebound from our optimization can be viewed as a worst-case estimate from a collection of models
that likely include the truth. In this sense, the bound is robust to model misspecification (to theextent of the neighborhood we impose). Our approach can also be viewed as a robustification of
traditional scenario analysis or stress testing. Rather than testing the impact of risk metric due tospecific shocks (by running the model at different scenarios or parameter values, say), our boundis capable of capturing the impact to the risk metric due to any perturbation of the model in themuch bigger non-parametric space. We will illustrate these interpretations in several examples suchas calculating loss probabilities and conditional value-at-risk.
The approach and the formulations that we discuss in this paper have roots in areas such as
economics, operations research and statistics. The notion of optimization over models appears
in the context of decision-making under ambiguity, i.e., when parts of the underlying model is3
uncertain, the decision-maker resorts to optimizing decisions over the worst-case scenario, resulting
in a minimax problem where the maximization is over the set of uncertain models. In the caseof probabilistic models, this machinery involves an optimization over distributions, which is theframework we study in this paper. In stochastic control, such an approach has been used inderiving best control policies under ambiguous transition distributions [54, 52, 37]. In economics,
the work of two Nobel laureates L. P. Hansen and T. J. Sargent [35] studies optimal decision
making and its macroeconomic implications under ambiguity. Similar ideas have been used inquantitative finance including portfolio optimization [29, 30], as well as in physics and biology
applications [3, 4]. In operations research, optimization over probabilistic models dates back to[60] in the context of inventory management, and has been used in [8, 9, 62] when some moment
information is assumed known. The literature of distributionally robust optimization, which has
been growingly active in recent years, e.g., [21, 32, 66, 38, 50], studies reformulations and efficientalgorithms to handle uncertainty in probabilistic assumptions in various stochastic optimization
settings. The specific type of constraints we study, namely a neighborhood defined via statisticaldistance, has been used in [6, 48, 42, 43, 36] under the umbrella of so-called φ-divergence, which inparticular includes the Kullback-Leibler divergence or the relative entropy that we employ heavily.
Other distances of recent interest include, e.g, the Wasserstein distance [28, 15]. In the context
of extreme risks relevant to actuarial applications, [14, 5, 10, 23] study the use of distances, e.g.,
Renyi divergence, to capture model uncertainty in the tail, and [44, 49] study the incorporationof tail shape information, along with the shape-constrained distributional optimization literature[56, 34, 47, 63]. In the multivariate setting, [26, 58, 64, 24, 57, 27] study analytical solution andcomputation methods for bounding worst-case copulas. [20] investigates the use of moments andentropy maximization in portfolio selection to hedge downside risks of a mortality portfolio. Lastly,
the calibration of the neighborhood size in distributionally robust optimization has been investigated
via the viewpoints of, e.g., hypothesis testing (e.g., [7]), empirical likelihood and profile inference(e.g., [65, 46, 11, 25, 41, 33]) and Bayesian (e.g, [7]). For consistency, throughout this paper,
we will adopt the terminology from the operations research literature and call our methodology
distributionally robust optimization.
Our contribution in this paper is to bring the combination of the ideas in the above areas

4
to the attention of the actuarial community. In addition, given that risk assessment is of special
importance for actuaries, we also discuss the implications of the distributionally robust methodology
in the setting of tail events beyond these past literatures. We will illustrate our modeling approachand interpretations, how it can be used in actuarial problems, and solution methods employing a
mix of analytical tools, simulation and convex optimization. We choose to present examples that are
relatively simple for pedagogical purposes, but our goal is to convince the reader that the proposed
methodology is substantially general.

2

Basic Distributionally Robust Problem Formulation

Suppose that in a generic application, one is interested in computing a performance measure in theform of an expected value of a function of underlying risk factors, i.e., Etrue (h (X)), where X is arandom variable (or random vector) taking values in Rd , and h : Rd → R is a performance function.
The notation Etrue (·) denotes the expectation operator associated with an underlying “true” orcorrect probabilistic model, which is unknown to the actuary.
To quantify model error, we propose, in the most basic form, to use the following pair of
maximization and minimization problems, which we call the basic distributionally robust (BDR)
problem formulation:

min / max E (h (X))(1)

s.t. D (P ||P0 ) ≤ δ,

The first involves solvingwhere
D(P kP0 ) = EP


log

 Z


dP
dP
(X)
= log
(x) dP (x) ,
dP0
dP0

(2)

is the so-called Kullback-Leibler (KL) divergence [40] between P and P0 , with the integral in (2)taken over the region in which X takes its values and dP/dP0 is the likelihood ratio between theprobability models P and P0 . Optimization (1) has a decision variable P , and the expectation in
its objective function is with respect to P . Finally, δ should be suitably calibrated (chosen as small

5
as possible) to guarantee that
D(Ptrue kP0 ) ≤ δ.

(3)

The KL divergence defined in (2) plays an important role in information theory and statistics (withconnections to concepts such as entropy [18] and Fisher information [19]; it is also known as therelative entropy whose properties have been substantially studied. The KL divergence measuresthe discrepancy between two probability distributions, in the sense that D(P kP0 ) = 0 if and only
if P is identical to P0 (also see discussion below). Thus, the constraint in (1) can be viewed as aneighborhood around P0 in the space of models, where the size of the neighborhood is measured by
the KL divergence.
The reason for selecting δ depicted above is the following. First, by choosing δ so that inequality(3) holds, we guarantee that Ptrue is in the feasible region of (1). This translates, in turn, to thatthe interval obtained from the minimum and maximum values of (1) contains the true performancemeasure Etrue (h(X)). Second, δ is chosen as small as possible because then this obtained intervalis the shortest, signaling a higher accuracy of our estimate.
The probability space of X described above can be substantially more general, including infinite
dimensional objects such as when X denotes a stochastic process like Brownian motion. However,to avoid technicalities, we will confine ourselves in the framework X ∈ Rd . In fact, to facilitate
the discussion further, first we discuss the case where P is discrete finite. In this case, (1) can be
written as

min / maxMX

p (k) h (k)

(4)

k=1

s.t.
D(pkp0 ) =

M
X
k=1

M
X

p (k) = 1,

p (k) logp (k)p0 (k)

≤δ,

p (k) ≥ 0 for k ≥ 1 ,k=1

where the decision variable is the probability mass function {p(k) : k = 1, . . . , M } for some support
size M , and {p0 (k) : k = 1, . . . , M } is the baseline probability mass function.
Let us briefly discuss two key properties of the KL divergence that make formulation (4) ap-

6
pealing. First, by Jensen’s inequality, we have thatD(pkp0 ) = −MXp (k) log

k=1

p0 (k)p (k)≥ − logMX
k=1

p0 (k)
p (k)
p (k)

!
= 0,

and D(pkp0 ) = 0 if and only if p (k) = p0 (k) for all k. In other words, the D (·) allows us tocompare the discrepancies between any two models, and the models agree if and only if there is nodiscrepancy. That is, in principle, the space of decision variables can include any distributions in
the form {p (k) : 1 ≤ k ≤ M }, without confining to any parametric form. Thus the framework is
capable to include model misspecification in a non-parametric manner.
It should be noted that D (·) is not a distance in the mathematical sense, because it does notobey the triangle inequality. However, and as the second key property of the KL divergence, D(pkp0 )is a convex function of p (·). To see this, first note that the function l (x) = x log (x) is convex on(0, ∞). Second, observe that if {p (k) : 0 ≤ k ≤ M } and {q (k) : 0 ≤ k ≤ M } are two probability
distributions and α ∈ (0, 1), then, because of Jensen’s inequality, for any fixed k,lαp (k) + (1 − α) q (k)p0 (k)




≤ αl

p (k)
p0 (k)

+ (1 − α) lq (k)p0 (k)

.

Thus, we conclude
M
X




αp (k) + (1 − α) q (k)D(αp + (1 − α) qkp0 ) =p0 (k) lp0 (k)
k=1




M 
X
p (k)
q (k)
≤
+ (1 − α) p0 (k)lαp0 (k)lp0 (k)
p0 (k)
k=1= αD(pkp0 ) + (1 − α) D(qkp0 ).

Consequently, BDR problem (4) is a convex optimization problem with a linear objective function.
These types of problems have been well-studied and are computationally tractable using operations
research tools.
Later in the paper we will discuss the implications of KL divergence as a notion of discrepancy
and propose other notions and constraints. For now, we argue that there are direct variations of the
BDR problem formulation that can be easily accommodated within the same convex optimization

7
framework. For this discussion, it helps to think of p(k) as the mortality distribution at each age (sok = 1, . . . , 100 for instance). Suppose that an actuary is less confident in the estimate for ptrue (k)for small values of k; that is, assume that p0 (k) ≈ ptrue (k) for large values of k, but the actuary is
uncertain about how similar p0 (k) is to ptrue (k) for small values of k (this arises if the mortality
estimate is more credible for some age groups for instance). Then we can replace the first inequality
constraint in (4) by introducing a weighting function {w (k) : 1 ≤ k ≤ M } with w (k) > 0 which isincreasing, thus obtainingMX
w (k) p (k) log

k=1

p (k)
p0 (k)

≤δ.

(5)To understand why w (·) should be increasing in order to account for model errors from misspecifying
p0 (k) for small values of k, consider the following exemplary case. Suppose for some k0 > 0, we
have 1) w (k) = ε > 0 for k ≤ k0 and ε small, and 2) w (k) = 1 for k > k0 . Then observe thatthe constraint (5) is relatively insensitive to the value of p (k) for k ≤ k0 . Therefore, the convex
optimization program will have more freedom for the p(k) on k ≤ k0 to jitter and improve theobjective function without having a significant impact on feasibility.
Another variation of the BDR formulation includes moment constraints. For instance, suppose
additional information is known for the expected time-until-death for individuals who have a particular underlying medical condition and are at least 30 years old (from a series of medical studies,say). Using such information one might impose a constraint of the form E(X|X ≥ 30) ∈ [a− , a+ ]for a specified range [a− , a+ ], or equivalently E(XI(X ≥ 30)) ∈ [P (X ≥ 30)a− , P (X ≥ 30)a+ ](throughout the paper, we use I (A) to denote the indicator of A, i.e., I (A) = 1 if A occurs and= 0 if not). Using this information, one can add the constraints

a−

M
X
k=30

p (k) ≤MX

p (k) k ≤ a+

k=30

M
X

p (k) ,

(6)k=30

which are linear inequalities, and therefore the resulting optimization problem is still convex andtractable. In Section 9 we will continue discussing other constraints that can inform the BDR
formulation with alternate forms of expert knowledge.
At this point, several questions might be in order: How do we solve the BDR optimization
problem and its variations? How can we understand such solution intuitively? What is the role8of the constraint in (1)? How do we choose δ? How do we extend the methodology to deal withpossibly multidimensional distributions? Our goal is to address these questions throughout the restof this paper.

3

Solving the BDR Formulation

This section describes how to solve (4) and then transitions toward the more general problem (1)in an intuitive way. We concentrate on the problem of maximization; the minimization counterpart
is analogous, and we will summarize the differences at the end of our discussion.

3.1

The Maximization Form

We introduce Lagrange multipliers to solve the convex optimization problem (4). The Lagrangiantakes the formg (p (1) , ..., p (M ) , λ1 , λ2 ) =

M
X
k=1

p (k) h (k) − λ1

M
X
k=1


p (k) logp (k)p0 (k)

!


−δ

− λ2

M
X

!
p (k) − 1 .

k=1

The Karush-Kuhn-Tucker (KKT) [16] conditions in our (convex optimization) setting characterize
an optimal solution. Denoting an optimal solution for (4) as {p+ (k) : 1 ≤ k ≤ M }, and the corre+
sponding Lagrange multipliers as λ+1 and λ2 , the KKT conditions are as follows (we use “+” todenote an optimal solution for the maximization formulation, as opposed to “−” for the minimiza-9tion counterpart, which we shall discuss momentarily as well):




∂g
p+ (k)+ ++p+ (1) , ..., p+ (M ) , λ1 , λ2 = h (k) − λ1 1 + log− λ+
2 ≤ 0,∂p (k)
p0 (k)
k = 1, . . . , M,(7)
p+ (k)∂g++p+ (1) , ..., p+ (M ) , λ+1 + log
1 , λ2 = p+ (k) h (k) − λ1∂p (k)

p+ (k)
p0 (k)

− λ+
2 = 0,

k = 1, . . . , M,


M
M
X
X
p+ (k)
p+ (k) log≤ δ,
p+ (k) = 1, p+ (k) ≥ 0 for 1 ≤ k ≤ M,
p0 (k)(8)+
λ+1 ≥ 0, λ2 is free,

(10)

k=1

λ+
1

(9)

k=1

M
X
k=1


p+ (k) log

p+ (k)
p0 (k)

!


−δ

= 0.

(11)

The relations (7) and (8) correspond to the so-called stationarity conditions, (9) the primal feasibility, (10) the dual feasibility, and (11) the complementary slackness condition.
DefineM+ = argmax{h (k) : 1 ≤ k ≤ M },i.e., M+ is the index set on which h (·) achieves its maximum value. Also, denote h∗ = max{h(k) :1 ≤ k ≤ M } as the maximum value of h(·). We analyze the solution by dividing it into two cases,1depending on whether log P0 (X∈M
≤ δ:
+)


1
Case 1: log P0 (X∈M
≤ δ.
+)Considerp+ (k) =p0 (k) I (k ∈ M+ )P
= P0 (X = k | X ∈ M+ ) .
j∈M+ p0 (j)(12)

for k ∈ M+ , and p+ (k) = 0 otherwise, i.e., p+ (·) is the conditional distribution of X given X ∈ M+ .
Note that this choice of p+ (·) gives
M
X
k=1


p+ (k) log

p+ (k)
p0 (k)




= log

1
P0 (X ∈ M+ )


≤δ

+
∗
so that (9) holds. Moreover, choosing λ+
1 = 0 and λ2 = h , we have (7), (8), (10) and (11) all hold.

Thus, (12) is an optimal solution in this case.

10Case 2: log

1
P0 (X∈M+ )



> δ.

Consider settingh (k) − λ+
1




p+ (k)
1 + log− λ+
2 = 0, k = 1, . . . , M,
p0 (k)+++++
giving p+ (k) = p0 (k) exp((h(k) − λ+
1 − λ2 )/λ1 ). Denoting θ+ = 1/λ1 and ψ+ = λ2 /λ1 + 1, we

can write
p+ (k) = p0 (k) exp (θ+ h (k) − ψ+ ) .
To satisfyPM

k=1 p+ (k) = 1 in (9), we must have ψ+ = log

(13)

PM

k=1 p0 (k) exp (θ+ h (k)), i.e., ψ+ is thelogarithmic moment generating function of h(X) under p0 (·) at θ+ . Now, to satisfy (11), we enforce
θ+ to be the positive root of the equation
M
X
k=1


p+ (k) log

p+ (k)
p0 (k)

!
M
X
p
(k)
exp
(θ
h
(k))
h
(k)
0
+
k=1
− log
p0 (k) exp (θ+ h (k)) = δ,PM
k=1 p0 (k) exp (θ+ h (k))
k=1
(14)

PM


= θ+

where the first equality is obtained by plugging in the expression of p+ (·) in (13). We argue thatsuch a positive root must exist. Note that the left hand side of (14) is continuous and increasing(continuity is immediate, and monotonicity can be verified by somewhat tedious but elementary
differentiation). When θ+ → 0 in (14), the left hand side goes to 0. When θ+ → ∞, it becomesθ+ h∗ + O(exp(−cθ+ )) − θ+ h∗ + log

X

p0 (k) + O(exp(−cθ+ )

k∈M+= − logX

p0 (k) + O(exp(−cθ+ ))k∈M+

for some constant c > 0, by singling out the dominant exponential factor exp(θ+ h∗ ) and noticing
a relative exponential decay in the remaining terms in the expression in (14). But note thatP
− log k∈M+ p0 (k) + O(exp(−cθ+ )) > δ as θ+ → ∞ in our considered case. Thus we must havea positive root for (14). Consequently, one can verify straightforwardly that the choice of p+ (·) in(13) with θ+ > 0 solving (14) satisfies all of (7), (8), (9), (10) and (11).
We comment that Case 1 can be regarded as a degenerate case and rarely arises in practice,

11
while Case 2 is the more important case to consider. From (13) in Case 2, p+ (·) can be interpreted
as a member of a “natural exponential family”, also known as an “exponential tilting” distributionthat arises often in statistics, large deviations analysis and importance sampling in Monte Carlosimulation. On the other hand, the form of p+ (·) in (12) in Case 1 can be interpreted as a limit of(13) as θ+ → ∞, namely

p+ (k) = lim p0 (k) PM

j=1 p0 (j) exp (θh (j))

θ→∞

3.2

exp (θh (k))=

p0 (k) I (k ∈ M+ )P
= P0 (X = k | X ∈ M+ ) .
j∈M+ p0 (j)The Minimization Form

The minimization form of the BDR problem is analogous to the maximization one. In this case,the Lagrangian takes the formg (p (1) , ..., p (M ) , λ1 , λ2 ) =MXp (k) h (k) + λ1

k=1

M
X
k=1


p (k) logp (k)p0 (k)

!


−δ

+ λ2

M
X

!
p (k) − 1 .

k=1With the corresponding KKT conditions and similar discussion as before, we arrive at the two
analogous cases. Here, we define

M− = argmin{h (k) : 1 ≤ k ≤ M }.

Case 1: log

1
P0 (X∈M− )



≤ δ.

An optimal solution is given byp− (k) = P0 (X = k|X ∈ M− ) .

Case 2: log

1
P0 (X∈M− )



> δ.

An optimal solution is given byp− (k) = p0 (k) exp (−θ− h (k) − ψ (θ− )) ,ψ (θ− ) = log

M
X
k=1

p0 (k) exp (−θ− h (k)) ,12with θ− > 0 satisfying
PM
−θ−

3.3

k=1 p0 (k) exp (−θ− h (k)) h (k)− logPM
k=1 p0 (k) exp (−θ− h (k))

MX!
p0 (k) exp (−θ− h (k))

= δ.

k=1

Solutions for the General Formulation

The insights of the solutions obtained in Sections 3.1 and 3.2 can be extended to the general BDR
formulation (1). First, denoting E0 (·) as the expectation under the baseline model P0 (·), (1) canbe written equivalently as optimizing over all random variables Z (ω), where ω is an element of theunderlying outcome space Ω, in the form

max E0 (h (X) Z)

(15)

s.t.
E0 (Z log (Z)) ≤ δ
Z (ω) ≥ 0, for all ω, and EP0 (Z) = 1 .

The object Z(ω) is the likelihood ratio between P and P0 , which is well-defined since the KL
divergence is defined only for P absolutely continuous with respect to P0 . The solution to (15),denoted (Z+ (ω) : ω ∈ Ω), can be used to obtain the solution P+ (·) to (1), by defining P+ (·) via

P+ (X ∈ A) = E0 (Z+ I (X ∈ A)) .

Conceptually, the general problem formulation is not much different from the finite outcome case(4). It can be shown (e.g., [54, 17]) that exactly the same form of the optimal solution applies. For
instance, if f0 (·) is the density of X in Rd , thenZ+ =f+ (X)= exp (θ+ h (X) − ψ (θ+ )) ,f0 (X)

ψ (θ+ ) = log E0 (exp (θ+ h (X))) ,

13
with
θ+

E0 (exp (θ+ h (X)) h (X))− log (E0 (exp (θ+ h (X)))) = δE0 (exp (θ+ h (X)))

in case log (1/P0 (X ∈ M+ )) > δ, where
M+ = argmax{h (x) : x ∈ Rd }.

However, if log (1/P0 (X ∈ M+ )) ≤ δ, then the optimal solution is a probability distribution givenby
P+ (X ∈ A) = P0 (X ∈ A|X ∈ M+ ) .

3.4

Numerical Example of BDR Formulation

Let us assume that X is the time-until-death of an individual entering a whole life insurance contractthat pays $1 benefit to the beneficiaries at the time of death. Assume that the force of interest(continuously compounded) is given by r > 0. We are interested in the expected net present
value of the benefits paid, denoted as h (X) = exp (−rX). This expected value is thus given by
Etrue (h (X)) = Etrue (exp (−rX)). A more realistic example might consider stochastic interest
rates or a portfolio of different types of contracts, but here we stay with this simple example forclear illustration.
Consider the more concrete setting where X takes values only on the set {1, ..., 100}, the benefit
is paid in integer years, and the individual in consideration is 20 years old so that the individual
could die at ages 21, 22,..., 120. Now, imagine a situation in which we know that the individual
has a medical condition for which not much is known, so there might be substantial variability in
our estimation of the distribution of X. Or perhaps there is a significant likelihood that medical
advances might be expected to occur in the next 10 years or so, and therefore, if p0 (·) were calibrated
based on historical data, p0 (·) might simply be an incorrect model.
To illustrate numerically the BDR formulation (4), we set the baseline distribution p0 (k) fromthe recent static mortality table under the Internal Revenue Code ([2], Appendix, “unisex” column),
and consider the estimation of the worst-case expected payoff when the true mortality distribution
deviates from p0 (k) as described by (4).

14
Figure 1 shows the maximum and minimum values of (4) under different δ, when the force ofinterest is set at r = .015. The maximum value increases with δ, while the minimum value decreases
with δ, both at a decreasing rate. The effects on the maximum values seem to be stronger than theminimum values when the distribution changes from the baseline, as shown by a larger magnitudeof the upward movement as δ increases. Note that the bounds, especially the upper one, grow quite
wide as δ increases. This indicates that under no other information, the worst-case impact of model
uncertainty can be quite significant. This impact can be reduced if one adds auxiliary information,either from additional data source or from expert knowledge (e.g., moment equalities or inequalitiesin (6) and in the subsequent Section 9).
Figures 2, 3 and 4 show the shapes of the distributions giving rise to the maximum and minimum
expected payoffs, compared with the baseline distribution, at δ = 1, 2 and 5 respectively. We cansee a sharp concentration of mass close to 0 for the maximal distribution, as higher chance of sooner
death leads to a larger expected present value of payoff. In contrast, more mass is located toward
older ages for the minimal distribution, leading to a reduction in the expected present value of
payoff. As δ increases, the accumulation of the mass at 0 for the maximal distribution tends to beheavier, as well as the shift of mass to older ages for the minimal distribution. The mass shift forthe maximal distribution seems to be more dramatic than for the minimal one, which is consistent
with Figure 1 that shows a wider bound for the maximization problem.
δ=1

Worst/Best−case Expected Present Value of Payoff

0.06
0.9

0.8

0.04

Optimal(Best−case)
0.6

Optimal(Worst−case)Original

Probability Mass

Expected Present Value

0.7

Distribution
P0
Pk−
Pk+

0.5
0.02

0.4

0.3

0.00

0

1

2

3

4

5

6

7

δ

Figure 1: Robust estimates against δ

0

25

50

75

100

k

Figure 2: Optimal probability mass function at δ = 1

15
δ=2

δ=5
0.15

0.06

0.10

Probability Mass

Probability Mass

0.04

Distribution
P0
Pk−
Pk+

Distribution
P0
Pk−
Pk+

0.05
0.02

0.00

0.00

0

25

50

75

100

0

25

k

50

75

100

k

Figure 3: Optimal probability mass function at δ = 2Figure 4: Optimal probability mass function at δ = 5

4Distributionally Robust Analysis of Dependence

We consider a variation of the BDR formulation that is suitable for quantifying the impact of
dependence. We revisit the example of whole life insurance in Section 3.4 as an illustration, but
now we assume that a couple is interested in a contract that pays $1 benefit at the time of the first
death between the couple. In this case, the payoff equals h (X, Y ) = exp (−r min(X, Y )), where X
and Y are the time-until-death of the individual and his or her spouse, respectively. We assume
that both spouses are 20 years old at the time of signing the risk premium valuation, so that thepair (X, Y ) is supported on the set {1, ..., 100} × {1, ..., 100}.
We are interested in the potential impact on the estimated actuarial net present value of thebenefits due to unaccounted dependence structures. Suppose that, this time, enough information
is available to estimate the distributions of X and Y marginally with relatively high accuracy, but
the joint distribution is difficult to estimate. That is, the marginals are known to be

Ptrue (Y = j) =100
X

P0 (X = i, Y = j) ,i=1

Ptrue (X = i) =100
X

P0 (X = i, Y = j) .

j=1

A joint baseline model p0 (i, j) = P (X = i, Y = j) can be conjectured (e.g., using an independent
marginal assumption).
To capture the above type of uncertainty, the distributionally robust (maximization) formulation

16
is given by

max

M X
M
X

p (i, j) h (i, j)(16)

i=1 j=1

s.t.
M X
M
X
i=1 j=1
MXp (i, j) log

p (i, j)p0 (i, j)≤ δ,p (i, j) = P0 (X = i) for all i,

j=1

M
X

p (i, j) = P0 (Y = j) for all j,

i=1

p (i, j) ≥ 0 for all i, j,

where we write (16) in a more general fashion that outputs the worst-case expected value E[h(X, Y )]over P (X, Y ) supported on {(i, j) : 1 ≤ i ≤ M, 1 ≤ j ≤ M }. The first constraint is a KL divergencefrom the baseline p0 (i, j) like in (1). The second constraint stipulates that the marginal distributions
are known to be P0 (X = i) and P0 (Y = j). In the joint whole life insurance example, M would be
set as 100.
We discuss how to solve (16), concentrating on the case analogous to Case 2 analyzed in the
BDR formulation in Section 3.1 (i.e., assuming the neighborhood size δ is not too big, which is themore interesting case).
It is worth mentioning that p = p0 is a feasible solution, and therefore we can verify the Slater
condition (see [16]). Such a condition guarantees that strong duality holds and that the KKT
conditions are necessary and sufficient for optimality. Consequently, commercial packages can beused to solve a convex optimization problem like (16) very quickly. Here, we will describe an
iterative procedure based on exponential tilting similar to that for problem (4). In particular, ananalogous argument in using the stationarity KKT condition in Case 2 there leads to an optimal
solution of the form

p+ (i, j) = P0 (X = i, Y = j) exp θ+ h (i, j) − αθ+ (i) − βθ+ (j) ,

17
where
X



exp θ+ h (i, j) − βθ+ (j) P0 (X = i, Y = j) = exp αθ+ (i) P0 (X = i) ,

(17)



exp θ+ h (i, j) − αθ+ (i) P0 (X = i, Y = j) = exp βθ+ (j) P0 (Y = j) .

(18)

j

X
i

where (17) and (18) guarantee the marginal distributions match the known values. The form above
motivates the following iterative scheme (see [22, 31]): Given a θ, defineexp (θh (i, j)) p0 (i, j).
pθ (i, j) = P Pkl exp (θh (k, l)) p0 (k, l)
Then iteratively update pθ (i, j) so that at each step, for each i, j = 1, . . . , M replace pθ (i, j) withpθ (i, j)P0 (X = i)Pl pθ (i, l)and then for each i, j = 1, . . . , M , replace pθ (i, j) withpθ (i, j)P0 (Y = j)Pk pθ (k, j)
until convergence. The iteration ensures the marginal distributions of the resulting probability massfunction are equal to P0 (X = i), i = 1, . . . , M , and P0 (Y = j), j = 1, . . . , M , respectively, leading
to a final pθ (·, ·). Then, we solve

min
θ

M X
M
X
i=1 j=1

pθ (i, j) h (i, j) −

1
θ

M X
M
X


pθ (i, j) log

i=1 j=1pθ (i, j)
p0 (i, j)




− δ

to get θ+ . An optimal mass function is then p+ (i, j) = pθ+ (i, j).

4.1

Numerical Example on Joint Mortality

We illustrate numerically the solution to optimization (16). Consider, as discussed above, theestimation of the expected payoff h(X, Y ) = exp(−r min(X, Y )) in the described setting. We set
the baseline joint distribution p0 (i, j) from the same mortality table used in Section 3.4, but now

18
assuming independence of the times-until-death of the couple, i.e. p0 (i, j) = P0 (X = i)P0 (Y = j).
This independent baseline model describes the simplest assumption made by insurers when pricing
life products, and we are interested in a robust estimate of the expected payoff when this assumption
is violated.
Figure 5 shows the worst scenario values of (16) under different values of δ, setting r = .015.
As depicted, there is a concavely increasing trend for the worst-case expected payoff as the level ofdependency deviates from independence. While the estimate is 0.0729 according to the independent
model, it can potentially rise to 0.0734 when the true model is misspecified from independence to alevel of δ = 0.002. Note that the magnitude of the changes in worst-case estimates is significantly
less dramatic than the example in Section 3.4. This is because here we have injected complete
marginal information, and the model is allowed to deviate only in terms of dependence structure.
From the small change magnitude, we can see that the marginal information is an important piecethat substantially reduces the freedom of deviation.
Figures 6 and 7 further show the approximate surfaces of the optimal joint probability massfunction at δ = 0.013 and the original independent joint probability mass function respectively.
Their differences are more clearly visualized in Figure 8. First, the maximum difference of the massfunctions is close to .01, which is quite small. This reconciles with the observation in Figure 5 thatthe deviations of the models are significantly confined due to known marginal information. Next,
in terms of the shape of the differences, the worst-case distribution falls below the baseline in aregion centered around (70, 70). This seems to be the region where the baseline probability massfunction peaks. The worst-case distribution spreads out this peak and allocates more masses toother regions to achieve a higher objective function. The positive difference seems clustered around
two regions and is interpreted as a consequence of the spreading period Note that the integrationof the differences must be zero, so both positive and negative differences must be present.

5

Distributionally Robust Rare-Event Analysis

Consider h (X) = I (X ∈ B) for some set B, where we recall I(·) as the indicator function. We areinterested in the probability Ptrue (X ∈ B) = Etrue [I(X ∈ B)], which is known to be small. We

19

Optimal Distribution
Worst−case Expected Present Value of Payoff
OptimalOriginal

●

●

●
●

0.15

0.0733●●0.10

0.0732

●

0.10
●

100
80

●

60

0.05

y

0.0731

Expected Present Value

0.15

●

0.0734●0.05

0.0730

40
0.00

20

0.0729

20

40

●

60

x

0.00000.00050.00100.0015

80

0.00

100

0.0020δ

Figure 6: Optimal probability mass function at δ =0.013

Figure 5: Robust estimate against δ

Difference between Optimal and Original Distributions
Initial Distribution
0.01

0.15

0.01
0.00

0.15

0.00
0.10

0.10

100

−0.01

100
80

y

60
−0.02

60

y

0.05

−0.01

80

40

0.05

−0.02

40
−0.03
0.00

20
20

40

60

x

80

100

20
20

60

x

0.00

Figure 7: Baseline probability mass function

40

80

100

Figure 8: Difference between optimal and baseline
probability mass functions

wish to solve the BDR maximization problem

max P (X ∈ B)

(19)

s.t. D (P ||P0 ) ≤ δ,

where a baseline distribution P0 of X has been suitably calibrated. Since we are in a rare-event
setting, it is reasonable to assume that log (1/P0 (X ∈ B)) > δ. Therefore, applying the general20
solution form of the BDR maximization problem in Section 3.3, a worst-case distribution is given
by P+ (·) defined such that for each set A,

P+ (X ∈ A) ==E0 (exp (θ+ I (X ∈ B)) I (X ∈ A))(20)
E0 (exp (θ+ I (X ∈ B)))exp (θ+ ) P0 (X ∈ A, X ∈ B)
P0 (X ∈ A, X ∈/ B)+,
exp (θ+ ) P0 (X ∈ B) + P0 (X ∈/ B) exp (θ+ ) P0 (X ∈ B) + P0 (X ∈/ B)

for some θ+ ∈ (0, ∞). To obtain a clearer interpretation of this worst-case distribution, let us write

P+ (X ∈ A) =P0 (X ∈ A)
(exp (θ+ ) − 1) P0 (X ∈ A, X ∈ B)+
exp (θ+ ) P0 (X ∈ B) + P0 (X ∈/ B) exp (θ+ ) P0 (X ∈ B) + P0 (X ∈/ B)

= α (θ+ ) P0 (X ∈ A|X ∈ B) + (1 − α (θ+ ))P0 (X ∈ A) ,where
α (θ+ ) =(exp (θ+ ) − 1) P0 (X ∈ B)
>0exp (θ+ ) P0 (X ∈ B) + P0 (X ∈/ B)and θ+ > 0 satisfiesθ+

exp (θ+ ) P0 (X ∈ B)
− log ((exp (θ+ ) − 1) P0 (X ∈ B) + 1) = δ.
(exp (θ+ ) − 1) P0 (X ∈ B) + 1(21)

Consequently, in simple words, in the context of distributionally robust performance analysisof rare-event probabilities, the worst-case measure is a (specific) mixture between the conditional
distribution of X given X ∈ B and the unconditional distribution, under the baseline measure(regardless of how it was calibrated). This allows the actuary to obtain bounds on rare-event
probabilities with the BDR formulation easily using only the baseline measure.
Moreover, note that the worst scenario value is given by

P+ (X ∈ B) =exp (θ+ ) P0 (X ∈ B)
exp (θ+ ) P0 (X ∈ B)=.
exp (θ+ ) P0 (X ∈ B) + P0 (X ∈/ B)
1 + (exp (θ+ ) − 1) P0 (X ∈ B)

(22)

Now, we investigate the asymptotic behavior of the worst scenario value (22) when we fix δ and
let P0 (X ∈ B) → 0. This asymptotics at a fixed model uncertainty level δ while the probabilityof interest is close to 0 is relevant from an applied standpoint in rare-event analysis. It models the21situation in which we have limited data, or postulate a specific imperfect parametric model P0 , yet
we want to estimate probabilities that are very small.
To perform the asymptotic analysis on P+ (X ∈ B), we first need to understand the behaviorof θ+ that satisfies (21). We carry out a heuristic development that follows and can be rigorized
by the argument in Theorem 1 in [10]. Suppose δ > 0 is fixed and P0 (X ∈ B) → 0. Let us
postulate that θ+ exp (θ+ ) P0 (X ∈ B) = η for some η > 0 that remains bounded away from zeroand infinity as P0 (X ∈ B) → 0. This implies that θ+ → ∞ as P0 (X ∈ B) → 0, and therefore
we have exp (θ+ ) P0 (X ∈ B) → 0 as P0 (X ∈ B) → 0. Consequently, (exp (θ+ ) − 1)P0 (X ∈ B) →0 and the left hand side of (21) is asymptotically equivalent to θ+ exp (θ+ ) P0 (X ∈ B), so thatθ+ exp (θ+ ) P0 (X ∈ B) ≈ δ as P0 (X ∈ B) → 0. This implies θ+ ≈ log(1/P0 (X ∈ B)). In turn,from (22), we get that
P+ (X ∈ B) ≈

δ
δ
.
≈
θ+
log (1/P0 (X ∈ B))

To understand decision-making implications, suppose that X models the potential losses of an
insurance company, and the baseline loss probability P0 (X > t) = exp (−Λ0 (t)) for some function
Λ0 (·), i.e., under P0 , X has a cumulative hazard function Λ0 (·). Now, say we wish to compute b
so that Ptrue (X > b) ≤ .005. Suppose that δ = .1, but the model Ptrue is not known. Then, basedon the worst-case distribution of our robust analysis above, we must choose b so that

P+ (X > b) ≈δ.1=
≤ .005,Λ0 (b)
Λ0 (b)(23)

which implies that b ≈ Λ−1
0 (20). For example, if X is exponentially distributed with unit meanunder P0 , so that Λ0 (b) = b, then this gives b ≈ 20. In contrast, if one trusts the exponential model
fully, then one would choose exp (−b) = .005, which yields b ≈ 5.3.
Paraphrasing, if b is the value of the statutory solvency capital needed to withstand losses withprobability at least .995, and if the actuary uses an exponential model with unit mean, then adeviation of .1 units measured in KL divergence between the assumed (exponential) model and theunknown reality might underestimate the statutory capital by a factor of about 20/5.3 ≈ 3.7.
Another important message is that the BDR formulation might provide very conservative estimates for rare-event probabilities when δ (i.e., the size of the uncertainty) is large relative to the22rare-event probability of interest. Suppose again that X is exponentially distributed with unit meanunder P0 , then (23) indicates that P+ (X > b) ≈ .1/b, giving a Pareto-type tail behavior under P+ .
In Section 9, we will discuss other types of distances that can give less conservative estimates forrare-event probabilities.

5.1

Numerical Example for Rare-Event Probability

Consider calculating the loss probability Ptrue (L > b) for a given reserve level b, where the loss model
is of the simple form L = X1 + · · · + Xn . Under the baseline model P0 , Xi ’s follow a multivariate
Gaussian distribution, i.e., (X1 , . . . , Xd ) ∼ N (µ, Σ) for some mean vector µ and covariance matrix
Σ, with specific distributional parameters given momentarily. Consequently,
P0 (L > b) = Φ

b − 10 µ
√
10 Σ1


,

where 1 is a vector of entries equal to one, 0 denotes transpose, and Φ (·) is the tail distributionfunction of a standard Gaussian random variable.
Now suppose that the actuary is uncertain about whether (X1 , ..., Xd ) follows a Gaussian model.
So we formulate the corresponding BDR problem (19) with the objective function being P (L > b).
We first check whether log (1/P0 (L > b)) ≤ δ, analogous to Case 1 in Section 3.1. If this is thecase, then the worst-case density of L is√φ (x − 10 µ) / 10 Σ1√
I(x > b),f+ (x) =P0 (L > b) 10 Σ1(24)

where φ (·) is the density of a standard Gaussian random variable, and P+ (L > b) = 1. To see this,note that h (l) = I (l > b), and (24) is precisely the conditional distribution of L given that L > b
under the model P0 , which is the solution depicted in Case 1 in Section 3.1.
Otherwise, if log (1/P0 (L > b)) > δ, then the worst-case density is given by

√
√
0
0
0

 exp(θ+ )φ((x−1 µ)/ 1 Σ1)/ 1 Σ1 for x > b
exp(θ+ )P0 (L>b)+P0 (L≤b)√√
f+ (x) =φ((x−10 µ)/ 10 Σ1)/ 10 Σ1
exp(θ+ )P0 (L>b)+P0 (L≤b) for x ≤ b

,

(25)

23
where θ+ > 0 is the root of the equationθ+exp (θ+ ) P0 (L > b)
− log(exp (θ+ ) P0 (L > b) + P0 (L ≤ b)) = δ.
exp (θ+ ) P0 (L > b) + P0 (L ≤ b)

These are obtained from (20) and (21), or analogously from Case 2 in Section 3.1. Alternately,using the interpretation of the associated exponential tilting discussed in Section 3.1, we have
 √x − 10 µ / 10 Σ1 exp (θ+ h (x))
 √∝ φ x − 10 µ / 10 Σ1 exp (θ+ ) I (x > b)
 √+φ x − 10 µ / 10 Σ1 I (x ≤ b) ,f+ (x) ∝ φwhich coincides with (25). Observe that the worst-case density is a weighted version of the original
one, thus largely preserving the shape in the corresponding regions. It assigns more (uniform)
weight to the region where x > b in a way that is as large as possible but preserves the validity of
the KL and the integrating-to-one constraints.
We illustrate the above numerically for n = 5, with means of Xi ’s being µ1 = 1.92, µ2 =
1.42, µ3 = 1.13, µ4 = 1.80, µ5 = 1.54 (five numbers generated uniformly between 1 and 2) andcovariance matrix


1.11−0.21 −0.42 −0.72

1.30







 −0.21 4.82
0.89 −0.31 −1.50 





Σ =  −0.42 0.89
1.05
0.61 −0.52 
,




 −0.72 −0.31 0.61
2.58 −0.45 
1.30 −1.50 −0.52 −0.45 1.91
which is randomly generated from a Wishart distribution with 5 degrees of freedom and an identity
scale matrix. Setting b = 10, the probability P0 (L > b) is now given by 0.23. Figure 9 shows theworst-case density under δ = 0.1 and the baseline density. Note that the worst-case density puts
more mass on the right side of b, and less on its left side, in order to boost the likelihood of a big loss.
In doing so, there is spike occurring exactly at b. This spike occurs because the adjustment leading
to the worst-case density attempts to maintain the shape of distribution as much as possible, so24that in each of the regions below and above b the density is multiplied by a constant weight. Finally,
Figure 10 shows, much like Figure 1, a concave growth pattern for the worst-case probability as δ
increases.

0.2
0.8

0.16

0.7

0.14

0.6

Worst-case probability

Density

0.18

0.12

0.1

baseline

0.08

worst-case

0.06
0.04
0.02

0.5
0.4
0.3
0.2

0.1

0
-5

0

5

10

15

20

0

0

Loss value

0.1

0.2

0.3

0.4

0.5

0.6

Delta

Figure 9: Approximate worst-case density under δ =0.1

6

Figure 10: Worst-case probability against δ

The Feasible Region in the BDR Formulation

This section discusses the selection of the parameter δ > 0. We present two main approaches: Oneis estimation using historical data. Another is to understand the choice of δ in terms of systematicstress testing. In the second approach, we will also establish the connection with conventional stresstesting approaches and parametric sensitivity analysis.

6.1

Data-Driven Estimation

For simplicity, we assume that the true model has a density ftrue (·), and the baseline model has adensity f0 (·). The most direct approach is to estimateZ
D (Ptrue ||P0 ) =log

ftrue (x)f0 (x)

ftrue (x) dx.

(26)

Some proposed procedures that guarantee convergence can be found in, e.g., [51, 55]. The
problem with this approach, however, is that the rate of convergence of the estimator depends on

25
the dimension of the underlying density. Intuitively, this is because this approach indirectly needs
to estimate density non-parametrically, which has a similar dependence on the dimension. An
alternate approach that is less conservative is based on empirical likelihood that we discuss next.
To explain, first suppose that we observe X1 , ..., Xn as an independent and identically distributed(i.i.d.) sample from X and form the empirical probability mass functionn1Xµn (x) =I (Xi = x) ,ni=1

which is uniform on the set {X1 , ..., Xn }. Now, for any given set of weights w = (w1 , ..., wn ) suchP
that wi ≥ 0 and ni=1 wi = 1, definevn (x, w) =n
X

wi I (Xi = x) .

i=1

In simple words, vn (·, w) is a probability mass function that assigns probability wi to the value Xi .
Then consider the so-called profile log-likelihood function [53]Rn (γ) = min D (vn (·, w) ||µn (·)) : Evn (·,w) (h (X)) =

Xx∈{X1 ,...,Xn }vn (x, w) h (x) =nXwi h (Xi ) = γ

i=1It turns out, under the null hypothesis that Etrue (h (X)) = γ and other mild conditions (including
V artrue (h (X)) < ∞), that we have
2nRn (γ) ⇒ χ21 .
That is, under the null hypothesis, 2nRn (γ) follows approximately a chi-square distribution with
one degree of freedom. This result is classical in the theory of empirical likelihood [53]. To make our
discussion as self-contained as possible, we provide formal derivation of the result in the appendix.
The connection with the BDR problem formulation can be established as follows. Assume P0
is built from the empirical distribution function of the observed data; that is, we use µn (·) as the
distribution of X under the model P0 . If we knew the specific value γ = Etrue (h (X)), it would
make sense to choose δ ≥ Rn (γ) to solve the BDR problem, because the worst-case distribution,





.

26
P+n , of the corresponding empirical BDR problem, namely

maxX

vn (x) h (x)(27)

x∈{X1 ,...,Xn }s.t. D (vn (·, w) ||µn (·)) ≤ δ ,will yield a worst scenario value, E+ (h (X)) =Pn+
k=1 wk h (Xk ) ≥ γ.

So, as n → ∞, the BDR

formulation will result in a correct upper bound for Etrue (h (X)).
Consequently, if one chooses P0 based on the empirical distribution, it is reasonable to select δ
so thatP (δ > Rn (γ)) ≈ P 2nδ > χ21 = .95,

(28)

thus providing a 95% confidence upper bound for γ. In other words, we can set δ = χ21,.95 /(2n),where χ21,.95 is the 95-percentile of χ21 . This approach hinges on the availability of a sufficiently
large sample size (in the same order as invoking the central limit theorem).

6.2

A Systematic Approach to Stress Testing

In the absence of enough data to perform a non-parametric calibration of P0 and δ as suggested
in the previous subsection, we suggest using the BDR formulation as an approach to performsystematic stress testing. The idea is to understand how the BDR optimization performs as wevary δ.
To motivate, consider the example of evaluating the expected shortfall or the conditional valueat-risk (C-VaR) of a given portfolio of risk exposures. C-VaR is the conditional expected size of the
deficit faced by the company, given the unlikely event of insufficiency of statutory solvency capital.
According to Solvency II, the capital should be sufficient to withstand losses with probability not
lower than .995 during a one-year time horizon [1]. To stress test C-VaR, it is customary toapply arbitrary shocks into the system (i.e., stressing the system by assuming that extreme eventsoccur). However, such shocks are typically selected in an arbitrary way. Moreover, in the presenceof multiple shocks affecting different risk factors, it might be difficult to argue that a particular
combination of shocks is more reasonable than an alternative combination of shocks.

27
In contrast, the approach that we study here can be used to stress test more systematically.
Suppose that an insurance company has a given baseline model, P0 , which has been calibrated using
a combination of past observations and expert knowledge. The model P0 is used to compute a givenrisk measure (or performance measure), say, C-VaR. Periodically, either by internal procedures ordue to regulatory constraints, the company is requested to perform stress testing. The result of
these stress tests is the determination of capital requirement, which might typically be higher thanthe C-VaR obtained under P0 .
Let us assume that such a stress-testing procedure has occurred multiple times in the company—say, n times—using the customary approach described before. So the company has built
certain experience of the increase in the capital requirement relative to the C-VaR determined by
the stress-testing procedure. Based on this experience, one could calibrate values δ1 , ..., δn so thatthe solution to the BDR formulation matches the increased capital requirement determined during the n sessions of stress testing. Ultimately, by choosing a single δ appropriately to embed allthese experiences, the BDR formulation automatically takes care of inducing the shocks that canpotentially cause the highest damage.
To demonstrate in a simple setting how the above could work, suppose that the shocks used instress-testing before are introduced by changing the parameters in P0 , which is assumed to have adensity lying in a parametric family, say fθ0 , where θ0 = (θ01 , . . . , θ0m ) denotes a set of m baseline
parameters. In the past stress-testing, the company changes the value of some θ0j to θ1j . Now, notethat the KL divergence between fθ and fθ0 , for any θ within a neighborhood of θ0 in the parametricspace, isZ
logfθ (x)fθ (x)dx = Eθ (log fθ (X) − log fθ0 (X))
fθ0 (x)where Eθ (·) denotes the expectation under fθ . By a Taylor series expansion, the above is approximately (see, e.g., [3, 4])
1Eθ log fθ (X) − log fθ (X) + ∇ log fθ (X)(θ0 − θ) + (θ0 − θ)0 ∇2 log fθ (X)(θ0 − θ)21=
(θ0 − θ)0 I(θ0 − θ),228
where ∇ log fθ and ∇2 log fθ denote the gradient and Hessian of log fθ , and I = −Eθ (∇2 log fθ (X))
is the Fisher information matrix. The last equality follows from the fact that the expected valueof the so-called score function ∇ log fθ , i.e, Eθ (∇ log fθ (X)), is 0. With this expression, suppose
that, in a simple setting, the company in a past period changes θ0j to θ1j with θ1j − θ0j = dj . Thismeans that by selecting δ = (1/2)Ij (dj )2 , where Ij is the j-th diagonal entry of I, the worst
scenario value of our BDR formulation gives a valid bound on the changed performance measureof interest from such a parametric stress test. Suppose that this has been done many times, sothat dj , j = 1, . . . , m, now denote the maximum perturbation difference used for each of the mparameters. Then choosing δ = maxj=1,...,m (1/2)Ij (dj )2 will give a robust stress test result thatincorporates all past experience.

7

Simulation-Based Solution Procedure

The optimization problems that we consider in the previous sections may sometimes be difficult
to solve, in the sense that an optimal model P+ is not expressible in closed form, especially in thecontext of general continuous distributions. For example, X may comprise a complicated functionor the convolution of many elementary variables. In this situation, one can resort to Monte Carlosimulation to approximate the optimization problem.
This approach is as follows. Assume that P0 is not analytically tractable or that we do not haveaccess to a closed-form expression of a density of X under P0 , but we can simulate i.i.d. sample
X1 , ..., Xn of X from P0 . Say we are interested in using BDR formulation (1) to bound Etrue (h (X)).
Using the re-expression (15), an empirical version of (15) takes the form given in (27), which we
write explicitly here as

max
s.t.

n
X
i=1
n
X
i=1

h (Xi ) wi
wi log (nwi ) ≤ δ ,

(29)
n
X

wi = 1, wi ≥ 0 for all 1 ≤ i ≤ n.

i=1

We then proceed to solve (29) using the tools for the discrete mass case in Section 3.1. Formulation(29) can be viewed as a so-called sample average approximation (SAA) (see [61, 39]), whose general29
idea is to replace expected value with empirical average in an optimization problem. Formulation(29) bears some differences with conventional SAA since its decision variable also changes with theMonte Carlo sample (depending on where the support of the sample is). [30, 31] use this approachto quantify model risk in some financial contexts.

7.1

Simulation-Based BDR Formulation: t-Copula Baseline Model

We illustrate the above simulation-based approach. Connecting with the example in Section 5.1,
suppose now that the loss is given by L = X1 + · · · + Xd , where (X1 , . . . , Xd ) has Gaussian marginal
distributions Xi ∼ N (µi , σi2 ) and the dependency is modeled by a t-copula. A t-copula is a multivariate distribution denoted by

t
−1
Cν,%
(u1 , . . . , ud ) = tν,% t−1
ν (u1 ), . . . , tν (ud ) ,(u1 , . . . , ud ) ∈ (0, 1)d ,where tν,% is the joint distribution function of a d-dimensional t-distribution with degree of freedom
ν, mean 0, and dispersion (or scale) matrix % = (ρij ) which is positive definite, and t−1
ν is thequantile function of a standard univariate t distribution with degree of freedom ν.
Then the distribution function of (X1 , . . . , Xd ) istf (x1 , . . . , xd ) = Cν,%(Φµ1 ,σ12 (x1 ), . . . , Φµd ,σ2 (xd )),dwhere Φµ,σ2 (·) denotes the distribution function of N (µ, σ 2 ).
It is difficult to evaluate P (L > b) in closed form, which motivates the use of Monte Carlosimulation. An unbiased estimate of P0 (L > b) can be obtained by outputting L = Φ−1(t (Z1 )) +µ ,σ 2 ν
1

1· · · + Φ−1(t (Zd )), where (Z1 , . . . , Zd ) is drawn from the multivariate t distribution tν,% . Repeatµd ,σd2 ν
this n times; say we get L1 , . . . , Ln .
To solve the empirical BDR formulation, we use the sampled L1 , . . . , Ln and combine with our
results in Section 5 as follows. We first check whether log(n/|{j : Lj > b}|) ≤ δ, where |{j : Lj > b}|

30is the cardinality of the set {j : Lj > b}. If this is the case, we let

for i such that Li > b 0

for i such that Li ≤ b

1|{j:Lj >b}|
+
wi =


.

This is the approximate worst-case distribution for L, which gives a corresponding approximate
worst scenario value 1.
Otherwise, if log(n/|{j : Lj > b}|) > δ, then output

wi =





exp(θ+ )exp(θ+ )|{i:Li >b}|+|{i:Li ≤b}|for i such that Li > b

1
exp(θ+ )|{i:Li >b}|+|{i:Li ≤b}|for i such that Li ≤ b,

where θ+ > 0 satisfies
θ+ exp (θ+ ) |{i : Li > b}|− log
exp (θ+ ) |{i : Li > b}| + |{i : Li ≤ b}|

1(exp (θ+ ) |{i : Li > b}| + |{i : Li ≤ b}|)
n


= δ.

The probability weights (wi )i=1,...,n on (Li )i=1,...,n form an approximation for a worst-case distribution for L, andexp (θ+ ) |{i : Li > b}|exp(θ+ )|{i : Li > b}| + |{i : Li ≤ b}|is the approximate worst scenario value P+ (L > b).
To illustrate numerically, we consider d = 5 and (X1 , . . . , X5 ) each following a Gaussian marginal
distribution with µ1 = 2.20, µ2 = 2.73, µ3 = 2.73, µ4 = 2.42, µ5 = 2.27 (five numbers generated
uniformly between 2 and 3) and σ1 = 0.92, σ2 = 0.39, σ3 = 0.11, σ4 = 0.56, σ5 = 0.33 (five numbers
generated uniformly between 0 and 1), respectively. We use a t-copula with a degree of freedom 10
and the following dispersion matrix:


8.19−0.92 −3.54

3.35

−3.96







 −0.92 6.09

1.07
1.37
−0.08





Σ =  −3.54 1.07
7.10
1.35 −1.70 
,




 3.351.37
1.35
3.14 −3.73 −3.96 −0.08 −1.70 −3.73 8.73

31
which is generated from a Wishart distribution with 5 degrees of freedom and an identity scale
matrix. We use n = 1, 000 to generate the baseline sample from the t-copula model. Setting b = 10,
the estimated loss probability is 0.232 with 95% confidence interval [0.206, 0.258]. The histogramsof the optimally weighted sample at δ = 0.1 and the baseline sample are plotted in Figures 11 and
12. As can be seen, more weights are put on the right side of b, in a uniform manner, to boost thelarge loss probability. We also plot the worst-case probability against δ in Figure 13, which, similar
to the example in Section 5.1, shows a concavely increasing pattern.

0.05
0.03

0.045
0.04

0.025

0.03

0.02

0.025

DensityDensity

0.035

0.02

0.015

0.015

0.01

0.01
0.005

0.005

11.75

0

0.85
1.25
1.65
2.05
2.45
2.85
3.25
3.65
4.05
4.45
4.85
5.25
5.65
6.05
6.45
6.85
7.25
7.65
8.05
8.45
8.85
9.25
9.65
10.05
10.45
10.85
11.25
11.65

11.35

10.95

10.55

9.75

10.15

9.35

8.95

8.55

8.15

7.75

7.35

6.95

6.55

6.15

5.75

5.35

4.95

4.55

4.15

3.75

3.35

2.95

2.55

2.15

0

Loss value

Loss value

Figure 11: Approximate worst-case density underδ = 0.1

8

Figure 12: Approximate original density

Robust Conditional Value at Risk

This section demonstrates how our distributionally robust analysis can be further developed tohandle risk-analytic problems involving optimization. For this discussion, we focus on the problem
of computing, say, the 95% C-VaR of a random variable L that represents potential losses in a given
year, i.e.,
C-VaRtrue (α) = Etrue (L − b|L > b) ,where Ptrue (L > b) = 1 − α with α = .95.

32

0.8

Worst-case probability

0.7
0.6
0.5
0.4

0.3
0.2
0.1
0
0

0.1

0.2

0.3

0.4

0.5

0.6

Delta

Figure 13: Worst-case probability against δ

It is well known [59] that if L has a continuous distribution under Ptrue , then

Etrue (max(L − θ, 0))
C-VaRtrue (α) = min θ +
θ
1−α


.

To robustify this calculation, we extend the BDR formulation discussed earlier in a natural way, tosolve an upper and a lower bound for C-VaRtrue (α) in the formE (max(L − θ, 0))
min / maxD(P ||P0 )≤δ min θ +
θ
1−α


,

(30)

where the outer optimization problems are taken over probability models P .
Problem (30) is challenging to solve in closed form, so it is natural to consider SAA along thediscussion in Section 7. To approximate
E0 (max(L − θ, 0))
C-VaR0 (α) = min θ +
θ
1−α

33
for a given baseline model P0 , we first simulate i.i.d. copies L1 , ..., Ln under P0 and then solven

X max(Li − θ, 0)\ 0 (α, n) = min θ + 1
C-VaR
θ
n
1−α

!
.

(31)

i=1

One property of this estimate is, from the theory of SAA, that under mild assumptions (for instance,
if L has a continuous density under P0 ),∗b (θ0 (n))\ 0 (α, n) ≈ C-VaR0 (α) + σ
C-VaR
Z,
n1/2

(32)

where Z is a standard Gaussian random variable,σb2(θ0∗ (n)) =

2
n 
1 X ∗
max(Li − θ0∗ (n) , 0)
\ 0 (α, n) ,
θ0 (n) +− C-VaR
n−1
1−α
j=1

and θ0∗ (n) is the solution obtained by solving (31). In simple words, σb2 (θ0∗ (n)) is the empiricalestimator of the variance of the objective function evaluated at the estimated optimal solution in(31).
Approximating (30) involves adding an outer optimization problem. That is, we still keep
L1 , ..., Ln simulated under P0 , but we consider

min / maxw1 ,...wn min θ +
θ

s.t.
n
X

n
X
i=1

max(Li − θ, 0)
wi
1−α

!
(33)

wi log (nwi ) ≤ δ

i=1
n
X

wi = 1, and wi ≥ 0 for all 1 ≤ i ≤ n.

i=1

The above discussion is similar to the BDR settings discussed before. However, there are also
some differences. One of these is that calibrating a suitable δ from the non-parametric, data-driven
view in Section 6.1 needs some modification in the C-VaR setting. It turns out that, with i.i.d. dataof size n, δ should be selected so that P 2nδ > χ22 = 1 − β, if we want the optimization output tobound the true C-VaR with (1 − β) confidence. Note that the degree of freedom of the chi-square

34
distribution has increased from one to two compared with the discussion in Section 6. The reason
for this change is the appearance of the inner minimization in problem (33), which introducesanother equality constraint in the empirical likelihood derivation outlined in the Appendix, in orderto capture the requirement that the derivative of the objective function with respect to θ should
vanish. This derivation is similar to that given in Section 6, and further details can be found in[45, 46]. Under additional assumptions, the calibration can be further tightened; see, e.g., [25].
Another point to note is that while the outer maximization formulation in (33) is a convex
program, the minimization one is not. Heuristic procedures, such as multi-start or alternating
minimization should be employed in the later case.

8.1

Numerical Example for C-VaR Estimation

We consider the problem of estimating C-VaR in which we assume that the loss L follows a standard
normal distribution under P0 . We set α = 0.9 and generate n = 1, 000 observations. Figure 14
shows the upper and lower robust bounds from (33) against different values of δ, showing how thewidth of the robust interval widens at a decreasing rate as δ increases.
Moreover, we also test the performance of the 95% (i.e., β = 0.05) confidence bounds using(33), by selecting δ such that P (2nδ > χ22 ) = 1 − β. We carry out the cases n = 50 and n = 100.
Table 1 reports the point estimate of the coverage probability, mean lower and upper bounds, andthe mean and standard deviation of the interval width for empirical likelihood, while Table 2 shows
the results using the classical SAA theory via (32). These quantities are obtained from repeating
the experiments multiple times, each time generating a new data set. The coverage probability,
for instance, outputs the proportion of repetitions in which the constructed interval covers thetruth. We see that empirical likelihood gives a higher and more accurate coverage, though the SAA
counterpart gives tighter and less varied interval width, at least for this example we conducted.
This demonstrates some benefits of using our BDR formulation with empirical likelihood calibration
over classical SAA, in generating intervals that cover the truth more frequently and hence are morerobust. On the downside, however, these more robust intervals pay the price of being generally
wider. The bottom line is that for a user who is risk-averse, the BDR approach with empirical
likelihood calibration should be preferable over SAA.

35

2.3371068
2.2871068
2.2371068

C-Var

2.1871068
2.1371068

upper

2.0871068

lower

2.0371068
1.9871068
1.9371068
1.8871068

0

1

2

3

4

5

6

7

8

9

Delta

Figure 14: Worst-case C-VaR against δ

n
50
100

Coverageprobability
0.90
0.94

Mean lower
bound
1.22
1.32

Mean upper
bound
2.33
2.26

Mean interval
width
1.11
0.94

Standard deviationof interval width
0.43
0.26

Table 1: Statistical performances of empirical likelihood for different sample sizes

n
50
100

Coverageprobability
0.86
0.84

Mean lower
bound
1.21
1.34

Mean upper
bound
2.26
2.05

Mean interval
width
1.05
0.71

Standard deviationof interval width
0.47
0.21

Table 2: Statistical performances of standard confidence interval of SAA for different sample sizes

9

Additional Considerations

We discuss two alternatives that can be used for robust performance analysis. The first one involves
the use of moment constraints, and the second one involving different notions of discrepancy.

9.1

Robust Performance Analysis via Moment Constraints

In some situations, it may not be possible to construct an explicit baseline distribution P0 . Alternatively, if information on moments is available, we might consider worst-case optimization under36such information, in the form
max E (h(X))
s.t.

E (vi (X)) ≤ αi , i = 1, . . . , s(34)

E(vi (X)) = αi , i = s + 1, . . . , m,where the maximization is over all probability models P (we focus on maximization here to avoid
redundancy). This is a general formulation that has m moment constraints, and vi (·) can represent
any function. For instance, for moment constraints involving means and variances, we can selectv1 (x) = x and v2 (x) = −x, v3 (x) = x2 , v4 (x) = −x2 , and α1 = µ, α2 = −µ, α3 = σ, α4 = −σ,and all constraints could be inequalities. There is a general procedure for solving these problems
which builds on linear programming. The most tricky part involves finding the support of the
distribution. Observe that, if the support of the worst-case distribution is known, then problem(34) is just a problem with a linear objective function and linear constraints, and is solvable using
standard routines.
Finding the support involves a sequential search. More precisely, the procedure for solving (34)is shown in Algorithm 1 (which borrows from, e.g., [9]).
We discuss Algorithm 1 in the following several aspects:
1. Interpretation: The output of the procedure is an exact optimal value of (34). The worst-caseprobability distribution is a finite-support discrete distribution on {x1 , . . . , xτ } with weights{pk1 , . . . , pkτ } obtained in the last iteration.
2. Comparison with the BDR formulation: Unlike the BDR formulation, (34) does not have abaseline input distribution to begin with.
3. Computational efficiency: Step 1 in each iteration of Algorithm 1 can be carried out bystandard linear programming solver, which can output both the optimal {pj } and the dualk }. Step 2 is a one-dimensional line search if X is one-dimensional.
multipliers {θk , π1k , . . . , πm

4. Minimization counterpart: For a minimization problem, simply replace h with −h in thePwhole procedure of Algorithm 1, except in the last step we output τj=1 h(xj )pkj .

37Algorithm 1 Generalized linear programming procedure for solving (34)
Initialization: An arbitrary probability distribution on the support {x1 , . . . , xL }, where l ≤ m+1,
that lies in the feasible region in (34). Set τ = l.
Procedure: For each iteration k = 1, 2, . . ., given {x1 , . . . , xτ }:1. Master problem solution: SolvePτh(xj )pjmax
Pj=1
τv (x )p ≤ αi , i = 1, . . . , ss.t.
Pτj=1 i j jv (x )p = αi , i = s + 1, . . . , m
Pτj=1 i j j
p
j=1 j = 1
pj ≥ 0, j = 1, . . . , τk } that satisfy
Let {pk1 , . . . , pkτ } be the optimal solution. Find the dual multipliers {θk , π1k , . . . , πmPkθk + Pm
i=1 πi vi (xj ) = h(xj ), if pj > 0, j = 1, . . . , τmθk + i=1 πik vi (xj ) ≥ h(xj ), if pj = 0, j = 1, . . . , τ
πik ≥ 0, i = 1, . . . , s

2. Subproblem solution: Find xτ +1 that maximizeskρ(x; θk , π1k , . . . , πm
) = h(x) − θk −mX

πik vi (x)i=1k ) > 0, then let τ = τ + 1; otherwise, stop the procedure, and {x , . . . , x }If ρ(xτ +1 ; θk , π1k , . . . , πm
1
τ
are the optimal support points, with {pk1 , . . . , pkτ } the associated weights.
After the last iteration, output
τXh(xj )pkj
j=1

38

9.2

Renyi Divergence as Discrepancy Notion

While we have presented our approach based on the KL divergence or relative entropy, other distancenotions can be used. For example, we could consider the BDR formulation using the so-called Renyi
divergence of degree α > 1, defined via
1
Dα (P ||P0 ) =log E0
α−1



dP
dP0

α 
.

As α → 1, we recover the KL divergence or relative entropy, i.e.,
 

dP
Dα (P ||P0 ) → D (P ||P0 ) = E log.
dP0
The corresponding distributionally robust formulation takes the form

min / max E (h (X)) s.t. Dα (P ||P0 ) ≤ δ,and the optimal solution is obtained roughly as follows. First, given θ1 , θ2 > 0, define
Z+ (θ1 , θ2 ) = max(θ1 + θ2 h (X) , 0)1/(1−α) ,with θ1 , θ2 chosen so that
E0 (Z+ (θ1 , θ2 )) = 1, E0 (Z+ (θ1 , θ2 )α ) = exp (δ (α − 1)) ,which comes from KKT conditions analogous to the one in Section 3.1. Then, a worst-case distribution is defined such that, for any set A,

P+ (X ∈ A) = E0 (Z+ (θ1 , θ2 ) I (X ∈ A)) .

See, e.g., [14, 5] for related derivations.
The use of Renyi divergence typically leads to less conservative estimates (as we discuss momentarily in the next subsection), but the parameter δ might be more difficult to estimate (e.g., the39empirical likelihood machinery discussed in Section 6.1 may not apply readily). Next, we consider
rare-event estimation to develop more intuition and contrast the effect of using Renyi divergencefrom KL.

9.2.1

Robust Rare-Event Analysis via Renyi Divergence

We consider the case in which h (X) = I (X ∈ B), so that we have
P+ (X ∈ B) = (θ1 + θ2 )1/(1−α) P0 (X ∈ B) ,andα/(1−α)E0 (Z+ (θ1 , θ2 )α ) = θ1

1/(1−α)E0 (Z+ (θ1 , θ2 )) = θ1

P0 (X ∈/ B) + (θ1 + θ2 )α/(1−α) P0 (X ∈ B) = exp (δ (α − 1)) ,P0 (X ∈/ B) + (θ1 + θ2 )1/(1−α) P0 (X ∈ B) = 1.

Letting (θ1 + θ2 )α/(1−α) = η/P0 (X ∈ B), and substituting in the previous display, we haveα/(1−α)θ11/(1−α)θ1(1 − P0 (X ∈ B)) + η = exp (δ (α − 1)) ,(1 − P0 (X ∈ B)) + η 1/α P0 (X ∈ B)1−1/α = 1.

As P0 (X ∈ B) → 0, since α > 1, we can select θ1 ≈ 1 and (θ1 + θ2 )1/(1−α) ≈ η 1/α /P0 (X ∈ B)1/α ,with η ≈ exp (δ (α − 1)) − 1, concluding that
P+ (X ∈ B) ≈ (exp (δ (α − 1)) − 1)1/α P0 (X ∈ B)1−1/α .

Let us revisit the example discussed at the end of Section 5. Assume that X is exponentiallydistributed with mean one under P0 . Select B = [b, ∞) and α = 2. Then we have
P+ (X > b) ≈ (exp (δ) − 1)1/2 exp (−b/2) .

In contrast, when using the KL divergence for the case of exponentially distributed X, we obtaina much lower rate of decay (of the form δ/b for b large; see (23)). In this sense, KL divergence40
induces a much more cautious robust methodology than Renyi divergence. To illustrate this lastpoint, we compare the BDR formulations using KL and Renyi divergences with α = 2. We use
h(X) = I(X > b) and an exponential baseline with rate 1, where we set b to be the 99-percentile ofthis baseline. We solve the BDR formulations using the technique in Section 5 and the discussion
above. Figure 15 shows that the bounds obtained from the Renyi divergence are tighter than KL in
this small probability estimation context. This is noticeable especially for the upper bounds where
the KL divergence gives values more than double in magnitude than Renyi as δ increases.

0.3

0.2

Probability

KL:Optimal(Best−case)
KL:Optimal(Worst−case)Original
RY:Optimal(Best−case)
RY:Optimal(Worst−case)

0.1

0.0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

δ

Figure 15: Robust estimates against δ using Renyi and KL divergences

10

Conclusions

We have discussed a systematic approach to quantify potential model errors based on the BDR
formulation, which is a convex optimization problem in the space of probability distributions. Thisapproach roots in areas including economics, statistics and operations research, and aims to computeworst-case bounds on the performance or risk measure when the model deviates from a baselinemodel in the KL sense. We have demonstrated several properties and advantages of this approach,including its non-parametric nature, generality in applying to various situations, and computational
tractability.
Building upon the BDR formulation, we have demonstrated how to extend our approach to41
applications that are relevant to actuarial practice, including studying the impact of dependence
in multivariate models, the impact of model uncertainty in rare-event estimation and in computing
risk-analytic measures that involve optimization such as C-VaR. We have also discussed methods
to calibrate the feasible region size in BDR formulations, from the viewpoint of empirical likelihoodand stress testing. Lastly, we have investigated additional ways to define the model uncertainty
region based on Renyi divergence or imposition of moment constraints.

11

Appendix: Approximating Distribution for the Size of the Feasible Region

Similar to the use of the KKT conditions for the solution of the BDR problem formulation, we
obtain that the optimal solution satisfies
exp (θh (Xi )),wi (θ) = Pn
j=1 exp (θh (Xj ))where θ satisfies

Pn

i=1 wi (θ) (h (Xi ) − γ) = 0, namely,Pn(θh (Xi )) (h (Xi ) − γ)i=1 exp
Pn= 0.
j=1 exp (θh (Xj ))(35)

Now suppose that indeed Etrue (h (X) − γ) = 0, and let us consider h̄ (x) = h (x) − γ. We canPrewrite (35)—after multiplying by nj=1 exp (θh (Xj ) − γθ) /n1/2 ) —as
n
1 X

n1/2


exp θh̄ (Xi ) h̄ (Xi ) = 0.

i=1Then let θ = η/n1/2 and perform a Taylor expansion to conclude that
n
1 X



h̄ (Xi )h̄ (Xi ) 1 + η 1/2 + · · ·
1/2
n
n
i=1


= 0.

(36)

42
Recall that the central limit theorem states the approximation in distributionn
1 X

n1/2 i=1h̄ (Xi ) ≈ σZ,

where σ 2 = V artrue (h (X)), and Z is standard Gaussian. Therefore, solving for η in (36) andignoring lower-order error terms, we obtain that

η≈−

1

Pn

n1/2 n1

Pi=1nh̄ (Xi )i=1 h̄ (Xi )Z

2 ≈ σ.

Now consider Rn (γ), which can be written as

Rn (γ)nX=
wi (θ) log
i=1
n
X

!

exp θh̄ (Xi )Pn−1 nj=1 exp θh̄ (Xj )!
n1X
exp(θh̄(Xi ))nwi (θ)h̄(Xi ) − log
i=1
i=1
1 Pn
θ n i=1 exp(θh̄(Xi ))h̄(Xi )
− log
1 Pn
i=1 exp(θ h̄(Xi ))
n

= θ

=
=

!
n
1X
exp(θh̄(Xi ))
n
i=1

P
P
θ n1 ni=1 h̄(Xi ) + nθ ni=1 h̄(Xi )2 + · · ·
P
θ2 Pn
2
1 + nθ ni=1 h̄(Xi ) + 2ni=1 h̄(Xi ) + · · ·


!2
n
n
n
2 X
2
X
X
θ
1
θ
θ
−
h̄(Xi ) +h̄(Xi )2 −h̄(Xi ) + · · · 
n
2n
2 n
i=1

n

≈

=

≈

i=1

n

θ2 X
θX
h̄(Xi ) +h̄(Xi )2
n
n
i=1
i=1

!2
n
n
n
2
X
X
X
θ
θ
1
1
− θ2
h̄(Xi ) −h̄(Xi ) −h̄(Xi )2 −
n
n
2 n
i=1
i=1
i=1


!
2
n
n
θ2  1 X
1X
2h̄(Xi ) −h̄(Xi ) 
2 n
n
i=1

≈

i=1

η2σ2
2n
Z2
,
2n

i=1

!2 
n
1X
h̄(Xi ) 
n
i=1

43
where the third last step follows by collecting terms up to θ2 .

Therefore, we conclude that

2nRn (γ) ⇒ χ21 , as indicated in our discussion leading to (28).

12

Acknowledgments(1) This project was supported by the Society of Actuaries (SOA) research grant entitled “Modelingand Analyzing Extreme Risks in Insurance.” The authors would like to thank the SOA for its
permission to publish some materials excerpted from our project reports [12, 13].
(2) We would like to thank the two referees for their careful reading and helpful feedback.

References[1] Directive 2009/138/ec of the European Parliament and of the Council of 25 November 2009 on
the taking-up and pursuit of the business of Insurance and Reinsurance (Solvency II). Available
at http://eur-lex.europa.eu/LexUriServ/LexUriServ.do?uri=OJ:L:2009:335:0001:0155:en:PDF .
[2] Updated static mortality tables for defined benefit pension plans for 2016.

Available athttps://www.irs.gov/pub/irs-drop/n-15-53.pdf .
[3] Arampatzis, G., M. A. Katsoulakis, and Y. Pantazis (2015a). Accelerated sensitivity analysis
in high-dimensional stochastic reaction networks. PloS one 10 (7), e0130825.
[4] Arampatzis, G., M. A. Katsoulakis, and Y. Pantazis (2015b). Pathwise sensitivity analysis intransient regimes. In Stochastic Equations for Complex Systems, pp. 105–124. Springer.
[5] Atar, R., K. Chowdhary, and P. Dupuis (2015). Robust bounds on risk-sensitive functionals via
rényi divergence. SIAM/ASA Journal on Uncertainty Quantification 3 (1), 18–33.
[6] Ben-Tal, A., D. Den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen (2013). Robust solutions of optimization problems affected by uncertain probabilities. Management Science 59 (2),341–357.
[7] Bertsimas, D., V. Gupta, and N. Kallus (2014). Robust SAA. arXiv preprint arXiv:1408.4445 .

44[8] Bertsimas, D. and I. Popescu (2005). Optimal inequalities in probability theory: A convex
optimization approach. SIAM Journal on Optimization 15 (3), 780–804.
[9] Birge, J. R. and J. H. Dulá (1991). Bounding separable recourse functions with limited distribution information. Annals of Operations Research 30 (1), 277–298.
[10] Blanchet, J., C. Dolan, and H. Lam (2014). Robust rare-event performance analysis withnatural non-convex constraints. In Proceedings of the 2014 Winter Simulation Conference, pp.
595–603. IEEE Press.
[11] Blanchet, J. and Y. Kang (2016). Sample out-of-sample inference based on wasserstein distance.
arXiv preprint arXiv:1605.01340 .
[12] Blanchet, J., H. Lam, Q. Tang, and Z. Yuan.
for actuarial applications.

Applied robust performance analysis

Society of Actuaries. Available at https://www.soa.org/research-

reports/2016/research-applied-robust-performance/ .
[13] Blanchet, J., H. Lam, Q. Tang, and Z. Yuan. Applied robust performance analysis for actuarial applications: A guide. Society of Actuaries. Available at https://www.soa.org/researchreports/2016/research-applied-robust-performance/ .
[14] Blanchet, J. and K. R. Murthy (2016a). On distributionally robust extreme value analysis.
arXiv preprint arXiv:1601.06858 .
[15] Blanchet, J. and K. R. Murthy (2016b). Quantifying distributional model risk via optimaltransport. arXiv preprint arXiv:1604.01446 .
[16] Boyd, S. and L. Vandenberghe (2004). Convex optimization. Cambridge University Press.
[17] Breuer, T. and I. Csiszár (2013). Measuring distribution model risk. Mathematical Finance.
[18] Cover, T. M. and J. A. Thomas (2012). Elements of information theory. John Wiley & Sons.
[19] Cox, D. R. and D. V. Hinkley (1979). Theoretical statistics. CRC Press.
[20] Cox, S. H., Y. Lin, R. Tian, and L. F. Zuluaga. Mortality portfolio risk management. Journalof Risk and Insurance 80 (4), 853–890.

45[21] Delage, E. and Y. Ye (2010). Distributionally robust optimization under moment uncertainty
with application to data-driven problems. Operations Research 58 (3), 595–612.
[22] Deming, W. E. and F. F. Stephan (1940). On a least squares adjustment of a sampled frequency
table when the expected marginal totals are known. The Annals of Mathematical Statistics 11 (4),427–444.
[23] Dey, S. and S. Juneja (2012). Incorporating fat tails in financial models using entropic divergence measures. arXiv preprint arXiv:1203.0643 .
[24] Dhara, A., B. Das, and K. Natarajan (2017). Worst-case expected shortfall with univariate
and bivariate marginals. arXiv preprint arXiv:1701.04167 .
[25] Duchi, J., P. Glynn, and H. Namkoong (2016). Statistics of robust optimization: A generalized
empirical likelihood approach. arXiv preprint arXiv:1610.03425 .
[26] Embrechts, P. and G. Puccetti (2006). Bounds for functions of multivariate risks. Journal ofMultivariate Analysis 97 (2), 526–547.
[27] Embrechts, P., G. Puccetti, and L. Rüschendorf (2013). Model uncertainty and var aggregation.
Journal of Banking & Finance 37 (8), 2750–2764.
[28] Esfahani, P. M. and D. Kuhn (2015). Data-driven distributionally robust optimization using
the Wasserstein metric: Performance guarantees and tractable reformulations. arXiv preprint
arXiv:1505.05116 .
[29] Glasserman, P. and X. Xu (2013). Robust portfolio control with stochastic factor dynamics.
Operations Research 61 (4), 874–893.
[30] Glasserman, P. and X. Xu (2014). Robust risk measurement and model risk. Quantitative
Finance 14 (1), 29–58.
[31] Glasserman, P. and L. Yang (2016). Bounding wrong-way risk in CVA calculation. Mathematical Finance, doi:10.1111/mafi.12141.

46[32] Goh, J. and M. Sim (2010). Distributionally robust optimization and its tractable approximations. Operations Research 58 (4-part-1), 902–917.
[33] Gotoh, J.-y., M. J. Kim, and A. Lim (2015). Robust empirical optimization is almost the sameas mean-variance optimization.
[34] Hanasusanto, G. A., V. Roitch, D. Kuhn, and W. Wiesemann (2015). Ambiguous joint chanceconstraints under mean and dispersion information. Technical report.
[35] Hansen, L. P. and T. J. Sargent (2008). Robustness. Princeton University Press.
[36] Hu, Z. and L. J. Hong (2013). Kullback-Leibler divergence constrained distributionally robust
optimization. Available at Optimization Online.
[37] Iyengar, G. N. (2005).

Robust dynamic programming.

Mathematics of Operations Re-

search 30 (2), 257–280.
[38] Jiang, R. and Y. Guan (2012). Data-driven chance constrained stochastic program. Mathematical Programming, 1–37.
[39] Kleywegt, A. J., A. Shapiro, and T. Homem-de Mello (2002). The sample average approximation method for stochastic discrete optimization. SIAM Journal on Optimization 12 (2), 479–502.
[40] Kullback, S. and R. A. Leibler (1951). On information and sufficiency. The Annals of Mathematical Statistics 22 (1), 79–86.
[41] Lam, H. (2016a). Recovering best statistical guarantees via the empirical divergence-based
distributionally robust optimization. Forthcoming in Operations Research. Available at arXiv
preprint arXiv:1605.09349 .
[42] Lam, H. (2016b). Robust sensitivity analysis for stochastic systems. Mathematics of Operations
Research 41 (4), 1248–1275.
[43] Lam, H. (2018). Sensitivity to serial dependency of input processes: A robust approach.
Management Science 64 (3), 1311–1327.

47[44] Lam, H. and C. Mottet (2017). Tail analysis without parametric models: A worst-case perspective. Operations Research 65 (6), 1696–1711.
[45] Lam, H. and E. Zhou (2015). Quantifying uncertainty in sample average approximation. In
Proceedings of the 2015 Winter Simulation Conference, pp. 3846–3857. IEEE Press.
[46] Lam, H. and E. Zhou (2017). The empirical likelihood approach to quantifying uncertainty insample average approximation. Operations Research Letters 45 (4), 301–307.
[47] Li, B., R. Jiang, and J. L. Mathieu (2016). Ambiguous risk constraints with moment andunimodality information. Available at Optimization Online.
[48] Love, D. and G. Bayraksan (2015). Phi-divergence constrained ambiguous stochastic programsfor data-driven optimization. Technical report, The Ohio State University, Columbus.
[49] Mottet, C. and H. Lam (2017). On optimization over tail distributions. arXiv preprint
arXiv:1711.00573 .
[50] Natarajan, K., D. Pachamanova, and M. Sim (2008). Incorporating asymmetric distributional
information in robust value-at-risk optimization. Management Science 54 (3), 573–585.
[51] Nguyen, X., M. J. Wainwright, and M. I. Jordan (2010). Estimating divergence functionals
and the likelihood ratio by convex risk minimization. Information Theory, IEEE Transactions
on 56 (11), 5847–5861.
[52] Nilim, A. and L. El Ghaoui (2005). Robust control of Markov decision processes with uncertaintransition matrices. Operations Research 53 (5), 780–798.
[53] Owen, A. B. (2001). Empirical likelihood. CRC press.
[54] Petersen, I. R., M. R. James, and P. Dupuis (2000). Minimax optimal control of stochastic
uncertain systems with relative entropy constraints. Automatic Control, IEEE Transactions
on 45 (3), 398–412.
[55] Póczos, B. and J. Schneider (2011). On the estimation of alpha-divergences.

48[56] Popescu, I. (2005). A semidefinite programming approach to optimal-moment bounds forconvex classes of distributions. Mathematics of Operations Research 30 (3), 632–657.
[57] Puccetti, G. and L. Rüschendorf (2012). Computation of sharp bounds on the distributionof a function of dependent risks. Journal of Computational and Applied Mathematics 236 (7),1833–1840.
[58] Puccetti, G. and L. Rüschendorf (2013). Sharp bounds for sums of dependent risks. Journalof Applied Probability 50 (01), 42–53.
[59] Rockafellar, R. T. and S. Uryasev (2000). Optimization of conditional value-at-risk. Journalof Risk 2, 21–42.
[60] Scarf, H., K. Arrow, and S. Karlin (1958). A min-max solution of an inventory problem. Studies
in the mathematical theory of inventory and production 10, 201–209.
[61] Shapiro, A., D. Dentcheva, and A. Ruszczynski (2014). Lectures on stochastic programming:
Modeling and Theory, Volume 16. SIAM.
[62] Smith, J. E. (1995). Generalized Chebychev inequalities: theory and applications in decision
analysis. Operations Research 43 (5), 807–825.
[63] Van Parys, B. P., P. J. Goulart, and M. Morari (2015). Distributionally robust expectation
inequalities for structured distributions. Optimization Online.
[64] Wang, B. and R. Wang (2011). The complete mixability and convex minimization problems
with monotone marginal densities. Journal of Multivariate Analysis 102 (10), 1344–1360.
[65] Wang, Z., P. W. Glynn, and Y. Ye (2015). Likelihood robust optimization for data-driven
problems. Computational Management Science 13 (2), 241–261.
[66] Wiesemann, W., D. Kuhn, and M. Sim (2014). Distributionally robust convex optimization.
Operations Research 62 (6), 1358–1376.