*Article* **Marriage and Individual Equity Release Contracts with Dread Disease Insurance as a Tool for Managing the Pensioners' Budget**

**Agnieszka Marciniuk \* and Beata Zmy´slona**

Department of Statistics, Faculty of Economics and Finance, Wroclaw University of Economics and Business, ul. Komandorska 118/120, 53-345 Wrocław, Poland; beata.zmyslona@ue.wroc.pl **\*** Correspondence: agnieszka.marciniuk@ue.wroc.pl

**Abstract:** In many countries around the world, demographic and civilization changes have brought about the phenomenon of aging societies. This phenomenon affects the economy, especially the pension and health care systems, causing difficulties in their financing. The implementation of a policy that would effectively manage the problem of the longevity risk is thus required. Using housing resources and private health insurance to improve retirees' living standards may serve this purpose. The instruments we propose comprise two variants of contracts: the first for a marriage, the second for an individual client. We analysed the cash flow in both the cases. The results suggest that the amount of cash flows related to reverse equity and dread disease insurance benefits depends on the spouse's economic status, age, and health conditions. The benefits of the two variants of the contract vary. This paper examines numerous strategies for selecting the type of the contract, taking into consideration the abovementioned factors.

**Keywords:** equity release; reverse annuity contract; critical health insurance; cash flow; financial protection of elderly

#### **1. Introduction**

Demographic changes have been caused by many factors. On the one hand, the developments and progress in medicine and the improved living conditions have extended life expectancy. On the other hand, the birth rate has been systematically decreasing. The consequence of these phenomena is the rapid increase in the over-65 population. In 2019, the percentage of people in this age group reached about 20.3% of the total European population (Eurostat 2020b). This group has become the fastest-growing segment of the population. It would not be an exaggeration to say that the world is turning grey.

The declining birth rate and increasing longevity worldwide have contributed to significant changes in the population's makeup and affected many branches of the economy, especially the labour and medical service markets. An increase in the elderly dependency ratio overburdens the work-age population. This indicates that in the future, pensions' funding gap will be a significant social problem. Moreover, the burden on health care triggered by the expenditures for the 65 and over population is becoming a large problem due to deteriorating health conditions, the presence of chronic diseases, and the need to provide long-term care for this social group. Ensuring financial security and health services to such a large elderly population is becoming a big challenge (Chen et al. 2021).

Aging implies changes in financial behaviour and the perception of many aspects of life, as well as the deterioration of health. Older adults' life satisfaction depends on such factors as their pensions, social support, and ability for self-care. Overall health and financial worries are significantly associated with life satisfaction (Borg et al. 2006, 2008). Therefore, it is also a crucial issue to provide financial peace and treatment options for the elderly (Korenman et al. 2021).

**Citation:** Marciniuk, Agnieszka, and Beata Zmy´slona. 2022. Marriage and Individual Equity Release Contracts with Dread Disease Insurance as a Tool for Managing the Pensioners' Budget. *Risks* 10: 140. https:// doi.org/10.3390/risks10070140

Academic Editors: Ermanno Pitacco, Annamaria Olivieri and Montserrat Guillén

Received: 24 February 2022 Accepted: 8 July 2022 Published: 12 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Many members of society, including retirees, own apartments or homes that they do not want to sell. On average, 70 % of people own their own properties in European Union (EU) countries (Eurostat 2020a); in Poland, this figure reaches almost 84 %, in Hungary almost 90 %, and in Romania over 95 %. Therefore, financial institutions in different countries have introduced equity release contracts for the retired (Hanewald et al. 2016), which provide an additional income in exchange for surrendering their real estates. Many variants of these contracts are available in European countries (also in Poland). The largest European market for equity release contracts is the United Kingdom (UK) market (Shao et al. 2015), while the largest world markets are the United States of America and Australia (Lee and Shi 2021).

Many studies and debates focus on equity release contracts. It is suggested that the housing equity of elderly homeowners may provide significant financial resources. The authors of the paper by (Mullings and Hamnett 1992) examined the participation of elderly persons in equity release schemes. In the article by (Toussaint and Elsinga 2009), attitudes toward the use of housing equity in several EU countries are examined. Based on 2005 data, researchers have observed that respondents used housing equity in various ways to help plan their finances. Since 2005, more people are aware that a possible method for funding the expenses associated with old age is through the use of equity release contracts (Marciniuk et al. 2020). Retirees are asset-rich and cash-poor; therefore, their housing wealth can be used to supplement their retirement income. Thus, the equity release contract is socially significant (IFoA 2019). The authors of the article (Sharma et al. 2020) point to various schemes for approaching equity release contracts, such as focussing on government subsidies for contracts. They conclude that the focus of government policy on equity release to tackle the challenges of an aging population is misplaced. However, treating equity release as an additional private source of income has many benefits.

In the richer countries of the world, an increase in demand for medical services has been observed. This increase is caused not only by demographic changes, such as the phenomenon of an aging population, but also by the occurrence of civilization diseases and technological progress in medicine. These additional factors result in increased health care expenses; consequently, it is impossible to provide high-level medical services without a mechanism of co-financing by patients. Two basic forms of co-payment can be considered: private payments, i.e., out of pocket expenditures, or private health insurance. Private health insurance can ensure access to additional funding, especially for critical or chronic illnesses. The funds obtained from this type of insurance can be used to improve the quality of life and treatment standards even in places where the state finances medical care. However, not everyone can afford to buy this insurance. Therefore, it has been necessary to introduce financial products combining different types of contracts (e.g., life insurance or equity release) with health insurance. Additional contracts are a source of financing the health insurance premium. The paper by (Webber 1993) presents an example of the financial contract, which meets the healthcare costs of the elderly. The authors emphasize that a major demand for the new products have been created by demographic changes, more widespread homeownership, and other political and cultural changes. A rapidly expanding market for products providing finance and healthcare for the elderly was observed, and its analysis indicates that the aged care insurance is both feasible and welfare-enhancing (Paolucci et al. 2015). Contracts that combine with critical health insurance are described in many papers. An example of a contract with life insurance is examined in the papers by (D˛ebicka and Zmy´slona 2018, 2019). The second example, considering a contract combined with a reverse annuity contract is described in article (D˛ebicka et al. 2015), which discusses an individual contract and one for a married couple (Zmy´slona and Marciniuk 2020). In the two mentioned articles, home equity release products combined with dread disease insurance are promoted as a potential solution to long-term care costs for the elderly.

The aim of this paper is to present two variants of these contracts, the first for both spouses and the second is a combination of two separate contracts for the wife and the husband. In Section 2, we describe the theoretical background of the contracts, and the multiplied Markov models are applied. Two varieties of contract with their cash flows are presented in Section 3. In Section 4, we analyse and compare the cash flows related to both options. The analysis results suggest the contract version affects the amount of cash flows connected to the reverse annuity contract and the dread disease insurance benefits, depending on the financial needs, health status, and spouses' age. The construction of the presented model is original. It is a continuation of the research conducted by the authors of the current work. Firstly, the individual model described in (D˛ebicka et al. 2015) was generalized to the marital model. In addition, the model was extended to include health benefit payments in two states (not only at the moment of diagnosis but also when the health condition worsens). Secondly, all analyses are performed based on the authors' own life expectancy tables for lung cancer patients in the critical stage described in (D˛ebicka and Zmy´slona 2016, 2019). Finally, all calculations are performed based on original computer programs in MATLAB.

#### **2. Contract Combining Reverse Annuity with Critical Health Insurance**

The new product consists of two elements, namely, the reverse annuity contract and the dread disease insurance. One of the equity release forms similar to the home reversion scheme is the reverse annuity contract (Marciniuk et al. 2020). The benefits derived from the new product can help improve the living conditions and the quality of life for the elderly as well as provide additional financial resources in case of a critical illness. The construction of this kind of contract is based on a model that describes an extended lifetime considering the risk of morbidity, severe disease, and life expectancy in the case of the critical stage.

Usually, both spouses are property owners and so the joint further lifetime of the wife and the husband should be considered. The spouses, who want to contract reverse annuity, receive a lifetime annuity in exchange for surrendering their real estate to a company. We consider two kinds of cash flows; the first group is connected with the reverse annuity contract. We distinguish a single premium and annuity benefits. This premium is not actually paid. It depends on the percentage *α* of real estate value. In the considered contract, the annuity benefits are paid with the option of selecting the last-surviving status (Marciniuk et al. 2020). This status means that the annuity benefits are paid every year while at least one of the spouses is alive and healthy. The annuity benefits depend on the pensioner's age and his/her further life expectancy and the value of the real estate W.

The second kind of cash flow concerns critical disease insurance. The survival time in case of falling ill with a severe disease depends on its course. We consider the health insurance premium that is paid every year that the insured party is healthy and alive. This insurance premium is paid as a percentage (1 − *β*) of the reverse annuity benefits. The health insurance benefits are firstly paid in case of a diagnosis. Moreover, the contract enables a benefit payment in the situation of deterioration of a patient's condition. This is important because the deterioration of health conditions in the case of a critical illness involves the necessity of care in the terminal state. The end-of-life heath expenditures are often huge. In many cases, families caring for terminally ill patients require financial support. For this reason, mild and critical health conditions are distinguished. In the critical stage the fatality rates grow rapidly, which decreases the probability of survival year by year. The analysis of the joint cash flows connected with the reverse annuity contract and health insurance requires the introduction of a multistate model, which is built separately for both spouses.

We have distinguished the following steps:

First step: Equity release—in exchange for surrendering their real estate, the clients at the retirement age receive a lifetime annuity (they do not pay a cash premium for that); the level of this theoretical annuity is computed in relation to the value of the house.

Second step: from the theoretical annuity, a portion 1 − *β* is paid to the insurer to finance health insurance; the remaining portion *β* is paid in cash to the annuitant as a lifetime annuity (level 1).

Third step: in case of a mild stage illness, the annuitant does not have to pay the health premium and therefore receives an annuity (level 2) equal to the initial theoretical annuity; at the moment of a diagnosis, additional benefits are paid through health insurance.

Fourth step: in case of the critical stage of an illness; an additional single benefit financed by health insurance is paid to the client (level 3).

#### *2.1. Multistate Markov Models*

The multistate model is described by the state space *S* = {1, 2, . . . , *N*} and a set of direct transitions between states T of the state space, where (*i*, *j*) denotes a direct transition from state *i* to *j* (*i* 6= *j*, *i*, *j* ∈ *S*). All the possible insured risk events up until the end of the contract can be described by this model (Pitacco 2014; Dhaene et al. 2017; Haberman and Pitacco 2018). *We consider the Markov model in which the stationary* Markov chain is described as a *discrete-time process* {*X*(*k*), *k* ≥ 0}*. The model describes the state and its changes for the contract at time k* = {0, 1, . . . , *n*}.

The construction of a contract combining the reverse annuity and dread disease insurance enables the introduction of two models. The first one is related to the reverse annuity contract and the latter with health insurance. The model connected with the reverse annuity contract considers only two states:


The model connected with health insurance considers the course of a dread disease. Thus, the following states are distinguished (Zmy´slona and Marciniuk 2020):


We assume that an insured could survive in the critical stages of illness during *h* years, and for that reason the fifth stage is extended by the introduction of states 5 (1) , 5 (2) , . . . ,5 (*h*) . The Markov model has a lack of memory; therefore, the presented model is extended by one additional state denoted by the number 3. This state describes a case when the insured fell mildly ill during the year. The multistate model related to the critical health insurance is presented in Figure 1 (Zmy´slona and Marciniuk 2020). The circles present the states, and the arcs describe the direct transitions between the states. *Risks* **2022**, *10*, x FOR PEER REVIEW 5 of 17

tion probabilities matrix for a person initially aged *x* (during the time interval [*x* + *k*, *x* + *k* 

[ ] ( ) 11 12

where <sup>11</sup> *q* stands for the probability of surviving denoted as *<sup>x</sup> <sup>k</sup> p* <sup>+</sup> and <sup>12</sup> *q* represents the probability of death denoted as *<sup>x</sup> <sup>k</sup> q* <sup>+</sup> . Considering the sex of the insured requires the introduction of two separate matrixes [ ] *<sup>x</sup>* **Q** and [ ] *<sup>y</sup>* **Q** , where *x* denotes the initial age

The transition probabilities matrix connected with the multistate model related to

0 0 0 0 ... 0 0 0 0 0 ... 0 0 0 0 0 0 ... 0 0 0 0 0 0 ... 0 ... ... ... ... ... ... ... ... ...

( ) () ( ) ()

*k q q*

0 0 0 0 0 0 ... 0 1 0 0 0 0 0 0 001

The transition probabilities matrixes in the female and male populations are denoted

The formulas for these elements are obtained based on the methodology of the multistate life table (increment-decrement table). The descriptions of the formulation and estimation processes are described in (Dębicka and Zmyślona 2016, 2019). These elements

(1)

11 12 14 16

*qq q q*

0 0 0 ... 0

23 24 26 33 34 36

*qq q qq q*

**Q** (2)

0 1 *<sup>x</sup> q q <sup>k</sup>* <sup>=</sup>

*ij q k P X xk jX xk i* = += = . We de-

*ij q k* . We define the following nonzero transi-

**Q** , (1)

12 1

55 56

45 46

*q q*

**Figure 1.** The multistate model for one spouse. depend on the following population rates: **Figure 1.** The multistate model for one spouse.

*2.2. The Probability Structure of the Model* 

note the elements of matrix **Q** ( ) *k* by ( ) [ ] *<sup>x</sup>*

of a woman and *y* the initial age of a man.

[ ]

by [*x*] **Q***h* and [ ] *<sup>y</sup>* **Q***<sup>h</sup>* , respectively.

*x h*

=

critical health insurance is given by

+ 1)) connected with the model for the reverse annuity contract

#### *2.2. The Probability Structure of the Model*

The probabilistic structure of the multistate model for the considered contract is given by the transition probability matrix *q* [*x*] *ij* (*k*) = *P*(*X*(*x*, *k* + 1) = *j*|*X*(*x*, *k*) = *i*). We denote the elements of matrix **Q**(*k*) by *q* [*x*] *ij* (*k*). We define the following nonzero transition probabilities matrix for a person initially aged *x* (during the time interval [*x* + *k*, *x* + *k* + 1)) connected with the model for the reverse annuity contract

$$\mathbf{Q}^{[\boldsymbol{x}]}(k) = \begin{pmatrix} q\_{11} & q\_{12} \\ 0 & 1 \end{pmatrix} \tag{1}$$

where *q*<sup>11</sup> stands for the probability of surviving denoted as *px*+*<sup>k</sup>* and *q*<sup>12</sup> represents the probability of death denoted as *qx*+*<sup>k</sup>* . Considering the sex of the insured requires the introduction of two separate matrixes **Q**[*x*] and **Q**[*y*] , where *x* denotes the initial age of a woman and *y* the initial age of a man.

The transition probabilities matrix connected with the multistate model related to critical health insurance is given by

$$\mathbf{Q}\_{h}^{[\boldsymbol{\chi}]}(k) = \begin{pmatrix} q\_{11} & q\_{12} & 0 & q\_{14} & 0 & 0 & \dots & 0 & q\_{16} \\ 0 & 0 & q\_{23} & q\_{24} & 0 & 0 & \dots & 0 & q\_{26} \\ 0 & 0 & q\_{33} & q\_{34} & 0 & 0 & \dots & 0 & q\_{36} \\ 0 & 0 & 0 & 0 & q\_{45^{(1)}} & 0 & \dots & 0 & q\_{46} \\ 0 & 0 & 0 & 0 & 0 & q\_{5^{(1)}5^{(2)}} & \dots & 0 & q\_{5^{(1)}6} \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \\ \\ 0 & 0 & 0 & 0 & 0 & 0 & \dots & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} . \tag{2}$$

The transition probabilities matrixes in the female and male populations are denoted by **Q** [*x*] *h* and **Q** [*y*] *h* , respectively.

The formulas for these elements are obtained based on the methodology of the multistate life table (increment-decrement table). The descriptions of the formulation and estimation processes are described in (D˛ebicka and Zmy´slona 2016, 2019). These elements depend on the following population rates:

*qx*+*k*—the probability of death;

*ϕx*+*k*—the dread disease mortality rate;

*χx*+*k*—the dread disease incidence rate (the morbidity rate);

*ψx*+*k*—the percentage of patients diagnosed in the critical stage;

*ξx*+*k*—the probability of health deterioration to the critical state;

*d* (*i*,*j*) *<sup>x</sup>*+*k*—the fatality rate in the population of the critically ill.

The above first three values are calculated for the whole population and the remaining three for the patient population.

The formulas for the elements of the matrixes **Q***<sup>h</sup>* [*x*] are presented in Table 1 (D˛ebicka and Zmy´slona 2016, 2019). The matrix formulas for the male **Q***<sup>h</sup>* [*y*] population are obtained analogously to those for the female population.


**Table 1.** The formulas for the estimators of the transition probabilities for health insurance.

#### **3. Two Version of Combining Contract**

A property is usually owned by both spouses, so the contract should be concluded by both a husband and a wife. The spouses may sign either one marriage contract for which they will receive one annuity benefit, or two individual contracts, in which two annuity benefits are paid independently. In each of the abovementioned cases, health insurance benefits are paid independently for the husband and the wife. Both solutions have advantages and disadvantages and provide a varied number of benefits. The choice of the type of contract will have an impact on the formation and differentiation of the reverse annuity and health insurance benefits.

#### *3.1. Marriage Contract*

The spouses sign one contract in which, in exchange for surrendering their rights to the property, they receive an annuity for life. A percentage *α* of the real value of the property *W* funds the benefit. Then, they pay two separate premiums for private health insurance. The funds designated for paying health insurance form a part of the annuity premium. We define *β* ∈ [0, 1] as the reverse annuity parameter. This parameter describes the percentage of the paid annuity. The remaining part 1 − *β* describes the health insurance premium. The impact of this parameter on the benefits is described in (Zmy´slona and Marciniuk 2020). The annuity benefit is paid in advance for both spouses at the beginning of the year, if at least one of the spouses is alive.

The cash flows connected with the marriage version of the model are obtained in three steps. In the first step, the annuity benefits are calculated based on the multiple state Markov model, which describes the further lifetime of a couple. It consists of 4 states (1—both spouses are alive, 2—the wife is dead, 3—the husband is dead, and 4—both spouses are dead).

Let *pi*(*k*) = P(*X*(*k*) = *i*) and **P**(*k*) = (*p*1(*k*), *p*2(*k*), *p*3(*k*), *p*4(*k*)) *T* be the probabilities of the remaining process {*X*(*k*), *k* ≥ 0} at states *k*, for *k* = 1, 2, 3 . . .. The matrixes describing the process duration in each state are shown in Table 2 for the wife and the husband.

**Table 2.** The process duration matrixes for a marriage reverse annuity contract.


The transition probability matrix for the 4-state Markov model (for a marriage) is defined as the following Hadamard product

$$\mathbf{Q}^{X\circ Y}(k) = \mathbf{H}^X(k) \circ \mathbf{H}^Y(k),\tag{3}$$

where *x* denotes the age of the wife and *y* the age of the husband. An element of matrix (2) equals the following product *q* (*X*◦*Y*) *ij* = *h X ij* · *h Y ij* for each *i* and *j*.

The details connected with this model are presented in (D˛ebicka and Marciniuk 2014). The formulas are determined based on the matrix notations described in (D˛ebicka 2013).

The course of the process is described by the matrix D. This matrix is defined as the following formula

$$\mathbf{D} = \begin{pmatrix} \mathbf{P}^T(0) \\ \mathbf{P}^T(1) \\ \cdots \\ \mathbf{P}^T(n) \end{pmatrix} \in \mathbb{R}^{(n+1)\times(N\cdot N)} \,\prime \tag{4}$$

where **P**(0) = (1, 0, 0, . . . , 0) *<sup>T</sup>* <sup>∈</sup> *<sup>R</sup>* (*N*·*N*) denotes the initial distribution vector, **P** *T* (*t*) = **P** *T* (0) *t*−1 ∏ *k*=0 **Q**(*X*◦*Y*) (*t*), where *N* is the state number and *n* is a contract duration.

This course takes into consideration both spouses' survival times.

The annuity benefit is determined by means of the following formula (D˛ebicka et al. 2022)

$$\ddot{b} = \frac{a\mathcal{W}}{\mathbf{M}^T \left(\mathbf{I} - \mathbf{I}\_{n+1}\mathbf{I}\_{n+1}^T - \mathbf{I}\_1\mathbf{I}\_1^T\right)\mathbf{D}(\mathbf{S} - \mathbf{J}\_4) + 1},\tag{5}$$

where **M***<sup>T</sup>* = 1, *v*, *v* 2 , . . . , *v n* ∈ *R <sup>n</sup>*+<sup>1</sup> denotes a discounting factor, which relates to the interest rate. The auxiliary vectors are used, where **I***<sup>t</sup>* = 0, 0, . . . , 1 *t*+1 , 0, . . . , 0*<sup>T</sup>* ∈ *R* (*n*+1)×1 , **S** = (1, 1, . . . ., 1) ∈ *R* (*N*·*N*)×1 , **J***<sup>i</sup>* = 0, 0, . . . , 1 *i* , 0, . . . , 0*<sup>T</sup>* ∈ *R* (*N*·*N*)×<sup>1</sup> and **I** is an identity

matrix with n+1 rows and columns. In this case, *N* = 2. The formula for the reverse annuity contract benefit in the classical notation is presented in (e.g., Marciniuk et al. 2020).

In the second step, the health insurance benefits are estimated. The health insurance premium is separately payable for the husband and wife. The further lifetime of a couple is described by the model presented in Figure 1. The part of the pension allocated to the payment of health insurance is divided between the spouses proportionally to the *R* parameter, which denotes a fraction of the men in the population who fall ill with the disease that is subject to dread disease insurance within a year. This parameter determines the risk of morbidity connected with a critical disease in the population of men. The decomposition of the annuity benefit is represented by the formula .. *b* = *β* .. *b* + *R*(1 − *β*) .. *b* + (1 − *R*)(1 − *β*) .. *b*. Thus, the annual health insurance premium is given by the following formulas: for a husband *p<sup>Y</sup>* = *R*(1 − *β*) .. *b* and for a wife *p<sup>X</sup>* = (1 − *R*)(1 − *β*) .. *b*. The health premiums are paid separately by the spouses for when they are healthy. Health insurance benefits are paid individually to the spouses depending on their health condition (at the moment of the diagnosis of a critical illness and in the event of deterioration of health), for a wife *c<sup>X</sup>* and for a husband *cY*. The formula for calculating the value of the health benefits is given by (compare Appendix A)

$$c = \frac{\mathbf{M}^T \mathbf{D} \mathbf{i} \mathbf{a} \mathbf{g} (\mathbf{C}\_{\text{out}} \mathbf{D}^T) \mathbf{S}}{\mathbf{M}^T \mathbf{D} \mathbf{i} \mathbf{a} \mathbf{g} \left(\mathbf{C}\_{\text{in}}^{(1)} \mathbf{D}^T\right) \mathbf{S}}. \tag{6}$$

Matrix **D,** describing the probability structure, is calculated on the basis of the matrix of transition probabilities given by (2) for the husband and wife separately. Thus, the number of columns in matrix **D** is reduced to *N*. The matrix of all premiums paid during the term of the contract is denoted as **C***out* ∈ *R* (*n*+1)×*N*. The matrix **C** (1) **in** ∈ *R* (*n*+1)×*<sup>N</sup>* equals 1 in the second and the fourth columns.

In summary, during the duration of the contract, spouses receive the following benefits:


Therefore, the spouses receive annuity benefits reduced by the health insurance premium (if they are healthy). In the event of a diagnosis or a deterioration of their health state(s), the spouses receive health insurance benefits individually. From the moment the diagnosis is made, the spouse does not pay the health insurance premium.

#### *3.2. Two Individual Contracts*

An individual version of the contract implies that a married couple signs two separate contracts. In this case, the premiums and benefits connected with health insurance and the reverse annuity contract are paid separately. The benefits of the reverse annuity contract are financed and determined based on the part of the property value owned by a given spouse (usually 50%) and is equal to *πII* = (0.5 · *W*) · *α*. The annuity benefits are calculated separately for a husband .. *b <sup>Y</sup>* and for a wife .. *b <sup>X</sup>*. The sum of the annuity benefits in an individual version is higher than the annuity benefit paid under the marriage contract. However, in the individual case, the annuity benefits are paid independently to each of the spouses only until their death. Thus, in the event of the death of one of the spouses, the surviving partner will be deprived of the part of the dead spouse.

Within individual contracts, the spouses may separately determine a part of an annuity benefit, which will be spent on the health insurance premium by setting the value of the reverse annuity parameter *β*. Therefore, the amounts of the health insurance premium are equal to *p II <sup>Y</sup>* = (1 − *βY*) · .. *b <sup>Y</sup>* and *p II <sup>X</sup>* = (1 − *βX*) · .. *b <sup>X</sup>* for a husband and a wife, respectively. In the next step, the health insurance benefits are estimated for a husband *c II Y* and a wife *c II X* .

The cash flows related to the individual version of the contract are separately modelled for a husband and a wife. The probability structure describing the further lifetime of a wife is denoted by **Q***<sup>h</sup>* [*x*] , whereas for a husband it is denoted by **Q***<sup>h</sup>* [*y*] (compare (2)).

The benefits are obtained separately for individual models for a wife and a husband in three steps. In the first step, the annuity benefits are obtained separately on the basis of the further lifetime for a wife and a husband using the following formula (compare D˛ebicka and Marciniuk 2014)

$$\ddot{b} = \frac{\alpha (0.5 \text{W})}{\mathbf{M}^T (\mathbf{I} - \mathbf{I}\_{n+1} \mathbf{I}\_{n+1}^T) \mathbf{D} \mathbf{J}\_1} \,\tag{7}$$

which is calculated separately for a husband and a wife. The matrix **D** is obtained on the basis of the transition probability matrix given by (1) separately for the spouses.

In the second step, the health insurance benefits are obtained using the model described in Figure 1. The health insurance benefits are calculated in a similar way as in the case of a marriage version of the contract using the Formula (6).

During two individual contracts, the spouses receive the following benefits:


The spouses receive annuity benefits reduced by the health insurance premium (if they are alive and healthy). The death of a spouse deprives the living partner of a part of their pension. Payments of dread disease insurance benefits are made in the event of the diagnosis of a serious illness or deterioration of the health state. From the moment the diagnosis is made, the spouse does not pay the health insurance premium.

#### **4. Results**

The retirees' family, economic, and health conditions determine their financial needs. The state of their health of has an important impact on their financial resources. Simultaneously, the elderly often have limited funds to buy an insurance policy against the risk of a dread disease. Therefore, the described contract may be an alternative to obtaining additional pensioner's income and funds for treatment and palliative care in the case of a dread disease.

A sample of numerical examples using two variants of the equity release contract is presented in this section. In the *marriage variant*, it is assumed that a married couple has a property valued at *W*. Together, they decide to sign the marriage reverse annuity and spend the percentage of 1 − *β* on withdrawing from both of their dread disease insurance plans. In the *individual variant*, the spouses also have a property valued at the same sum *W*. They decide to sign the individual reverse annuity contracts for 0.5*W* and allocate the percentage of 1 − *β* to buy dread disease insurance from their part of the benefit payment. In this case, they have two separate contracts.

The empirical examples are based on actual data. The examples concern the risk of lung cancer mortality and the morbidity of the Lower Silesia population (one of Poland's voivodeships). The calculations are made separately for male and female groups based on datasets provided by the National Cancer Registry for the Lower Silesia Region (Wojciechowska and Didikowska 2014) and the Lower Silesia Department of the National Health Fund (the public payer in Poland). The available, unpublished data covers the period between 2006 and 2011. The transition probabilities matrixes are obtained from the life expectancy tables described in the (D˛ebicka and Zmy´slona 2016). The maximum survival time in a critical condition was 4 years. A critical stage entails diagnosis with distant metastases or an inoperable life-threatening tumour (D˛ebicka and Zmy´slona 2019).

Let *x* denote a woman's age and *y* denote a man's age. The percentage of a real value of property *α* = 50 and the real value of a property is 100,000 euro. This assumption allows for the easy rescaling of the results obtained in case of different property prices. The discounting factor *v*, which is calculated from formula *v* = (1 + *i*) −1 , is closely connected with the analysed period of the disease and a constant long-term interest rate *i =* 0.0579. The interest rate *i* was estimated on the basis of actual Polish market data based on the yield to maturity on the fixed interest bonds and Treasury bills from 2008 (Zmy´slona and Marciniuk 2020). The percentage of men who fall ill with lung cancer in 2008, denoted by *R*, equals 0.691. We consider couples aged 65, 70, 75, 80, and 85, as well as mixed-age couples. We analyse a case when *β* = 0.99, which means that 1% of the annuity is intended for paying the premium. It is obvious that the levels of the premium may vary (Zmy´slona

and Marciniuk 2020). Annual marriage reverse annuity payments, single benefit payments of illness insurance, which could be paid twice, and premiums in euros are calculated using self-developed programs written in MATLAB. Table 3 presents all cash flows in both variants.


**Table 3.** Cash flows in two variants of equity release.

In both variants, an increase in all the benefits was observed. This growth is positively correlated with the spouses' age. Women pay lower premiums because the incidence of lung cancer is lower in the female population. Consequently, the critical illness insurance benefit is significantly higher for females than for the male population. It is obvious that in case of another disease the relation may be opposite. The percentage of the morbidity rate under consideration *R* determines the distribution of the premium. In the marriage variant, the total benefit payment of the equity release is lower than the same benefit in the individual variant. The differences in the benefit payments increase with the rise in the spouses' age and reach almost 30% for *x* = *y* = 85. The same relationship is observable for illness benefits; however, women pay an almost two-fold higher premium in the second option, thus receiving twice as much payment. This shows that women could pay half the premium for this particular disease, which means 0.5 · (1 − *β*), to receive the same chronic illness benefit payment. The critical illness insurance payment for men is comparable in both cases. This is due to the fact that the premium is divided according to the morbidity risk in the marriage contract. In both variants, the spouses receive a reverse annuity as a yearly payment, and they can allow themselves to buy critical illness insurance. The annuity income is not significantly reduced because the level of the annual premium is very low. Whereas the health insurance benefit allows for the acquirement of significant financial resources for health care, treatment, and improvement of the quality of life in the case of dread disease.

As shown by the above analysis, it can be assumed that the second variant is more favourable for the spouses. Therefore, let us analyse the relative increases in the number of total benefits payments received by spouses at 65, 75, and 85 years depending on their further life expectancy, which is presented in Figures 2–4.

*Risks* **2022**, *10*, x FOR PEER REVIEW 12 of 17

**Figure 2.** The relative increase in benefits in both variants for spouses aged 65. **Figure 2.** The relative increase in benefits in both variants for spouses aged 65. 14.73, 9.22, and 5.2 years, respectively.

**Figure 3. Figure 3.** The relative increase in benefits in The relative increase in benefits in both variants for spouses aged 75. both variants for spouses aged 75.

**Figure 4.** The relative increase in benefits in both variants for spouses aged 85. **Figure 4.** The relative increase in benefits in both variants for spouses aged 85.

It is visible that the similar age of the spouses' death determines the total benefit of the equity release contract, increasing up to 40%, and results in a higher dread disease insurance payment in the event of sickness. This means that the second variant of the

rectangle.

especially be observed for women, who become widows more often than men, and receive a higher marital annuity, even by over 40%. However, when the husband additionally suffers from a serious illness, in the considered case of lung cancer, the health insurance benefit is of a similar amount. The lower left rectangle shows situations when the spouses live longer than their future life expectancy. Then, the probability of receiving a higher benefit increases. When spouses are at 65, the higher payments of the second variant of the contract occur in 40.7% of cases, of which 61.3% concerns a situation when spouses live above their future life expectancy. These numbers increase with the rise in the spouses' age. The second variant's higher benefit occurs in 50.9% and 63.1% when spouses are at 75 and 85 years, respectively. Of these situations, 72.3% are placed in the left lower

The spouses are not very often at the same age; therefore, Tables 4 and 5 present the amount of the annuity and health insurance benefits in both variants for spouses of various ages. An exemplary five-year age difference between the spouses was assumed.

**Table 4.** Total benefits of equity release for various ages of spouses in two variants of contract.

**First Variant—Marriage Contract Second Variant—Individual Contracts**  *x X y* 65 70 75 80 85 *y* 65 70 75 80 85 65 4085 4394 4717 5006 5226 65 5074 5444 6039 6982 8460 70 4224 4629 5084 5526 5889 70 5545 5915 6510 7453 8931 75 4340 4839 5458 6115 6712 75 6240 6610 7205 8147 9625 80 4428 5008 5788 6718 7647 80 7272 7641 8237 9179 10,657 85 4493 5133 6050 7241 8612 85 8775 9145 9740 10,683 12,161

The dark colours indicate positive differences, i.e., a higher total benefit in the second variant. Lighter colours represent negative differences, which means a higher total payment in the first option. The graphs also show lines that illustrate the future life expectancy of people at the considered age. The future life expectancy of a 65, 75, and 85-year-old woman is 19.04, 16.54, and 6.05 years, respectively, and for men of the same age it is 14.73, 9.22, and 5.2 years, respectively.

It is visible that the similar age of the spouses' death determines the total benefit of the equity release contract, increasing up to 40%, and results in a higher dread disease insurance payment in the event of sickness. This means that the second variant of the contract is more favourable for spouses. When one of the spouses dies earlier than the other, and the second spouse lives a long time, the first variant is more beneficial. This can especially be observed for women, who become widows more often than men, and receive a higher marital annuity, even by over 40%. However, when the husband additionally suffers from a serious illness, in the considered case of lung cancer, the health insurance benefit is of a similar amount. The lower left rectangle shows situations when the spouses live longer than their future life expectancy. Then, the probability of receiving a higher benefit increases. When spouses are at 65, the higher payments of the second variant of the contract occur in 40.7% of cases, of which 61.3% concerns a situation when spouses live above their future life expectancy. These numbers increase with the rise in the spouses' age. The second variant's higher benefit occurs in 50.9% and 63.1% when spouses are at 75 and 85 years, respectively. Of these situations, 72.3% are placed in the left lower rectangle.

The spouses are not very often at the same age; therefore, Tables 4 and 5 present the amount of the annuity and health insurance benefits in both variants for spouses of various ages. An exemplary five-year age difference between the spouses was assumed.


**Table 4.** Total benefits of equity release for various ages of spouses in two variants of contract.

It is evident that the reverse annuity payments increase alongside the rise in the spouses' age in both cases. However, the benefits are higher in the second variant. The difference in both benefits increases for younger women and decreases for older women as the men's age increases. For example, when x = 65 and y = 85, the benefit is over 95% higher, and when men become older, it is only 41.2% higher, but it is still huge. It is not difficult to notice that the payment is lower for a younger wife and an older husband than contrariwise. In the case of men, age has a higher impact on the benefit. This relates to the fact that men live shorter lives than women.

The critical illness insurance payment for women is almost always twice as high in the second variant. Wives do not share the premium for this insurance with husbands. This benefit does not depend on the age of the second spouse in the second case. For a younger husband and an older wife, the dread disease insurance payment is higher in the first contract; thus, it is more favourable for men.


**Table 5.** Critical illness insurance payment depending on the age of spouses in two variants of contract.

#### **5. Discussion and Conclusions**

The article presents two variants of a contract that is a combination of an equity release and critical illness insurance. This new proposition could protect against the effects of longevity. The contractual cash flow and different scenarios were analysed. It cannot be clearly stated which of these contracts is better. Assuming a pessimistic variant, when both spouses fall ill and die quickly, neither option will be beneficial, although they will receive the health insurance benefit payment. If one spouse dies earlier and the other lives healthily longer than the life expectancy, then the marriage variant is more favourable. On the other hand, assuming an optimistic version, when both spouses live together in good health for a long time, particularly above the life expectancy, the individual variant is definitely a better choice. Of course, when concluding a contract, most people do not consider the risk of death. The knowledge of the genetic burden concerning a critical illness can help to more adequately adjust the client's contract. It was shown that both contracts offer high benefits payments, but for women, who become widows more often, the marriage contract is more beneficial. The calculations show a very significant gender impact on the amount of the benefits. In the considered examples, net premiums and net benefits were presented. All the contracts have additional costs. Combining arrangements into one lowers the costs; hence, the first variant seems to be a cheaper option than buying individual policies.

The introduction of critical illness insurance into the contract has caused two important constraints for the presented model. Firstly, the limited critical illness survival time (a fouryear period for lung cancer as described in the studied example) enforced the application of the periodic life tables. Secondly, since an estimation of the relationship between cancer risk for both the spouses seems impossible, we have assumed independence between the spouses' life expectancies. A justification of the possibility of such an assumption of independence is presented in (Zmy´slona and Marciniuk 2020).

The studied contracts could be used as part of longevity risk management policies as they provide an instrument of incentive for the retirees to use housing resources to improve their living conditions. In case of a dread disease, additional funds are likely to improve the quality of life during the treatment. The two variants of contracts, marriage and individual, could help pensioners to manage their home budget.

**Author Contributions:** Conceptualization, A.M. and B.Z.; methodology, A.M. and B.Z.; software, A.M.; validation A.M. and B.Z.; formal analysis, A.M. and B.Z.; investigation, A.M. and B.Z.; resources, B.Z.; data curation, A.M. and B.Z.; writing—original draft preparation, A.M. and B.Z.; writing review and editing, A.M. and B.Z.; visualization, A.M. and B.Z.; supervision, A.M. and B.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** The project was financed by the Ministry of Science and Higher Education in Poland under the programme "Regional Initiative of Excellence" 2019—2022 project number 015/RID/2018/19 total funding amount 10 721 040,00 PLN.

**Data Availability Statement:** Analysis are based on: 1. The Life Tables from D˛ebicka, J., Zmy´slona, B. (2016). Construction of Multi-State Life Tables for Critical Illness Insurance–Influence of Age and Sex on the Incidence of Health Inequalities. Sl ˛aski Przegl ˛ad Statystyczny, 14, 41–63. ´ https: //doi.org/10.15611/sps.2016.14.03; 2. The National Cancer Registry for the Lower Silesia Region from Wojciechowska, U., and Didkowska, J. (2014). Zachorowania i zgony na nowotwory zło´sliwe w Polsce. Krajowy Rejestr Nowotworów. Centrum Onkologii—Instytut im. Marii Skłodowskiej—Curie. www.onkologia.org.pl/raporty/. 3. The Life Table for the Lower Silesia from 2008 (followed from Statistical Offices.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

The formula for the value of the health benefits is given by

$$\mathcal{L} = \frac{\mathbf{M}^T \mathbf{D} \mathbf{i} \mathbf{a} \mathbf{g} \left( \mathbf{C}\_{\text{out}} \mathbf{D}^T \right) \mathbf{S}}{\mathbf{M}^T \mathbf{D} \mathbf{i} \mathbf{a} \mathbf{g} \left( \mathbf{C}\_{\text{in}}^{(1)} \mathbf{D}^T \right) \mathbf{S}} \tag{A1}$$

**Proof.** Let **C** = **C**in − **C**out, where **C**out ∈ *R* (*n*+1)×*<sup>N</sup>* denotes the matrix of all premiums paid during the term of the contract and **Cin** ∈ *R* (*n*+1)×*<sup>N</sup>* is the matrix of all the benefits, which are paid at state 2 and 4 in the amount *c*. Hence


We use the equivalence principle for the multistate insurance (cf. D˛ebicka 2013), which has the following form

$$\mathbf{M}^{\mathsf{T}} \mathbf{Diag} \left( \mathbf{C} \mathbf{D}^{\mathsf{T}} \right) \mathbf{S} = 0,$$

where **M***<sup>T</sup>* = 1, *v*, *v* 2 , . . . , *v n* ∈ *R n*+1 , **S** = (1, 1, . . . ., 1) ∈ *R* (*N*·*N*)×<sup>1</sup> and **<sup>D</sup>** <sup>∈</sup> *<sup>R</sup>* (*n*+1)×(*N*·*N*) describes the matrix probability (such as in Formula (4)).

Hence

$$\mathbf{M}^T \mathbf{D} \mathbf{i} \mathbf{g} \left( (\mathbf{C}\_{in} - \mathbf{C}\_{\text{out}}) \mathbf{D}^T \right) \mathbf{S} = \mathbf{0}$$

and

$$\mathbf{M}^{\mathsf{T}} \mathbf{Diag} \left( \mathbf{C}\_{\mathsf{out}} \mathbf{D}^{\mathsf{T}} \right) \mathbf{S} = \mathbf{M}^{\mathsf{T}} \mathbf{Diag} \left( \mathbf{C}\_{\mathsf{in}} \mathbf{D}^{\mathsf{T}} \right) \mathbf{S}$$

Due to **Cin** = *c* · **C** (1) **in** , we receive

$$\mathbf{M}^T \mathbf{Diag} \left( \mathbf{C}\_{\text{out}} \mathbf{D}^T \right) \mathbf{S} = c \cdot \mathbf{M}^T \mathbf{Diag} \left( \mathbf{C}\_{\text{in}}^{(1)} \mathbf{D}^T \right) \mathbf{S}$$

which completes the proof.

#### **References**


## *Article* **A Comparison of Macaulay Approximations**

**Stefanos C. Orfanos**

Department of Risk Management & Insurance, J. Mack Robinson College of Business, Georgia State University, Atlanta, GA 30302, USA; sorfanos@gsu.edu; Tel.: +1-404-413-7485

**Abstract:** We discuss several known formulas that use the Macaulay duration and convexity of commonly used cash flow streams to approximate their net present value, and compare them with a new approximation formula that involves hyperbolic functions. Our objective is to assess the reliability of each approximation formula under different scenarios. The results in this note should be of interest to actuarial candidates and educators as well as analysts working in all areas of actuarial practice.

**Keywords:** Macaulay duration; Macaulay convexity; net present value of cash flows

#### **1. Introduction**

Actuaries and actuarial science students at universities all over the world are familiar with approximation formulas for the present value of cash flow streams using some notion of cash flow duration or convexity. For example, the syllabus of Exam FM of the US-based Society of Actuaries includes the topic of approximations using the Macaulay and modified duration and convexity, while the UK-based Institute and Faculty of Actuaries in its material for exam CM1 mentions approximations derived from a Taylor series expansion.

Beside the academic and pedagogical interest in such approximation formulas, one may also consider the practical value in the management of interest rate risk. Although abundant computing power has enabled firms to implement elaborate immunization strategies that incorporate multi-factor stochastic interest rate models, non-parallel yield curve shifts, and complicated asset and liability characteristics, the restrictions posed by a simplistic valuation model are not unreasonable if rates remain historically low, yield curves stay relatively flat, and we can control the potential errors. Indeed, it may be helpful to know which approximation formula proves to be the most reliable, and to use it as a quick validation tool when time constraints preclude the use of a more sophisticated approach.

Alps (2017) describes a realistic scenario involving an investment actuary and her CEO, where the use of an approximation formula would be warranted or even necessitated. This is especially true in today's world of fast-changing rates, when companies have to react almost instantly to benchmark fund rates and quantitative tightening decisions by the Federal Reserve or other central banks.

In this note, we discuss several known formulas that use the Macaulay duration and convexity of commonly used cash flow streams to approximate their net present value, and compare them with a new approximation formula that involves hyperbolic functions. In addition to annuities, dividend stocks, and bonds, we also consider the cases of negative payments and embedded options to perform a deeper assessment. The notions of effective duration and convexity are defined in the next section and used to price the embedded options. Our objective is to measure the reliability of each approximation formula under different scenarios. As alluded to earlier, we only consider parallel interest rate shocks to flat yield curves.

**Citation:** Orfanos, Stefanos C. 2022. A Comparison of Macaulay Approximations. *Risks* 10: 153. https://doi.org/10.3390/ risks10080153

Academic Editors: Ermanno Pitacco and Annamaria Olivieri

Received: 13 June 2022 Accepted: 15 July 2022 Published: 29 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### *1.1. Literature Review*

The idea of using a bond's duration to approximate changes to its price goes back to Macaulay (1938). Some authors credit Fischer and Weil (1971) with the publication of the first duration–convexity approximation formula. Enhancements of that formula by controlling the missing higher-order terms and incorporating passage of time were announced in Jarjir and Rakotondratsimba (2008, 2012), though the resulting formulas contain parameters that are unintuitive and hard to calibrate. A conceptually simpler formula was given in Barber (1995), and independently in Livingston and Zhou (2005) for the modified duration, which was subsequently generalized to a duration–convexity model in Tchuindjo (2008). Further work in Barber and Dandapani (2017) considered negative-yielding bonds, and Johansson (2012) added passage of time. A separate duration– convexity formula appears in Alps (2017) and is applied to an empirical study of basic immunization strategies in Nie et al. (2021), while a very recent paper by Barber (2022) further generalizes a duration–convexity approximation by introducing an additional 'compounding' parameter. Finally, traditional approximations have been implemented in statistical analysis packages; see Lee (2021) for R code.

#### *1.2. Notation*

We denote the net present value of a cash flow stream by *P*. The interest rate *r* is annualized and continuously compounded (i.e., force of interest). ∆*r* is the change in interest rates from the initial value *r*<sup>0</sup> to *r*. Finally, the annual discount factor *v* is by definition equal to *e* −*r* .

Throughout the remainder of this paper and for convenience, assume *r*<sup>0</sup> = 1.6%, which is approximately the yield on the 10-year T-bond at the beginning of this year.

#### **2. Materials and Methods**

Recall that the *Macaulay duration* of a stream of cash flows {*CFt<sup>j</sup>* } *n j*=1 being paid at future times {*tj*} *n j*=1 is defined by

$$d = \frac{\sum\_{j=1}^{n} \mathbb{C} F\_{t\_j} \ v^{t\_j} \ t\_j}{\sum\_{j=1}^{n} \mathbb{C} F\_{t\_j} \ v^{t\_j}} = -\frac{dP/dr}{P} \prime$$

while its *Macaulay convexity* is given by

$$\mathcal{c} = \frac{\sum\_{j=1}^{n} \mathcal{C} F\_{t\_j} \, v^{t\_j} \, t\_j^2}{\sum\_{j=1}^{n} \mathcal{C} F\_{t\_j} \, v^{t\_j}} = \frac{d^2 P / dr^2}{P} \,.$$

We do not consider the modified duration here because the Macaulay duration has a more intuitive interpretation (being the 'average' timing of the cash flows) and tends to result in tighter approximations for non-negative rates.

In the case of bonds with embedded options, it will be necessary to price the value of the option using a simple Black model. Recall that the pricing formula for, say, a European call option with expiration at time *t* and strike *K* is

$$V = v^{\mathfrak{t}}(P\Phi(d\_1) - K\Phi(d\_2))$$

where Φ represents the standard normal c.d.f. and the quantities *d*1,2 are given by

$$d\_{1,2} = \frac{\ln(P/K)}{\sigma\sqrt{t}} \pm \frac{\sigma\sqrt{t}}{2}$$

with *σ* the bond price volatility. We also need more flexible measures of bond duration and convexity. To that end, define the *effective duration* by means of

$$d\_{\varepsilon} = -\frac{P(r\_0 + \Delta r) - P(r\_0 - \Delta r)}{2P\_0 \,\Delta r}$$

and the *effective convexity* as

$$c\_{\varepsilon} = \frac{P(r\_0 + \Delta r) - 2P\_0 + P(r\_0 - \Delta r)}{P\_0 \ (\Delta r)^2}.$$

#### *2.1. Fischer–Weil's Approximation*

This follows immediately from Calculus and the definitions above.

$$\frac{\Delta P}{P\_0} \approx -d\_0 \,\Delta r + \frac{c\_0}{2} \, (\Delta r)^2. \tag{1}$$

It is assumed that the Macaulay duration and convexity are computed at rate *r*0, hence the subscripts.

#### *2.2. Barber' 1995 Approximation*

Instead of the second-order Taylor polynomial of *P*, we consider the first-order Taylor polynomial of ln *P*, thus obtaining

$$\ln P \approx \ln P\_0 - d\_0 \,\Delta r\_\prime$$

thus

$$P \approx P\_0 \, e^{-d\_0 \, \Delta r} \,. \tag{2}$$

Unlike the first-order Taylor polynomial in *P* that has no convexity, the functional form of Barber's approximation bequeaths it with a certain degree of positive curvature. This leads to good approximation results whenever *c*<sup>0</sup> ≈ *d* 2 0 and poor performance for *c*<sup>0</sup> < 0.

#### *2.3. Tchuindjo' Approximation*

Similar to Barber's approximation, but involving the second-order Taylor polynomial of ln *P*

$$
\ln P \approx \ln P\_0 - d\_0 \,\Delta r + \frac{c\_0 - d\_0^2}{2} \, (\Delta r)^2 \,\text{\textdegree} \tag{3}
$$

from which one solves for *P*. The added quadratic term gives better results in cases where *c*<sup>0</sup> − *d* 2 0 is non-trivial, but may still introduce large errors whenever *c*<sup>0</sup> < 0.

#### *2.4. Alps' Approximation*

The approximation formula and its derivation can be found in Alps (2017). The central idea in the derivation of this approximation is to compute a Taylor polynomial for the current value of the cash flow stream at time *t* = *d*0. This choice results in high accuracy in situations where *d*<sup>0</sup> > 0, and less so for *d*<sup>0</sup> < 0.

We have rewritten it below in terms of continuously compounded interest rates.

$$P \approx P\_0 \, e^{-d\_0 \, \Delta r} \left( 1 + \frac{c\_0 - d\_0^2}{2} \left( e^{\Delta r} - 1 \right)^2 \right). \tag{4}$$

#### *2.5. Hyperbolic Approximation*

We have not encountered this approximation formula in the literature and we assume its derivation is presented here for the first time. Consider the homogeneous differential equation *P* <sup>00</sup> − *c*0*P* = 0 that mimics the definition of Macauley convexity given earlier in this note. Its general solution takes the form

$$P = a\,e^{\sqrt{c\_0}\,^r} + b\,e^{-\sqrt{c\_0}\,^r}$$

(do not worry for the time being about the case *c*<sup>0</sup> < 0.) Setting *P*(*r*0) = *P*<sup>0</sup> and *P* 0 (*r*0) = −*d*<sup>0</sup> *P*0, which is a reformulation of the definition of Macauley duration, one obtains the approximation

$$P \approx P\_0 \left( \frac{1}{2} \left( 1 - \frac{d\_0}{\sqrt{c\_0}} \right) e^{\sqrt{c\_0} \Delta r} + \frac{1}{2} \left( 1 + \frac{d\_0}{\sqrt{c\_0}} \right) e^{-\sqrt{c\_0} \Delta r} \right)^2$$

which can be rewritten as

$$P \approx P\_0 \left( \cosh(\sqrt{c\_0} \,\Delta r) - \frac{d\_0}{\sqrt{c\_0}} \sinh(\sqrt{c\_0} \,\Delta r) \right). \tag{5}$$

The well-known trig identities

$$\cosh(i\theta) = \frac{e^{i\theta} + e^{-i\theta}}{2} = \cos\theta \, \, \, \qquad \sinh(i\theta) = \frac{e^{i\theta} - e^{-i\theta}}{2} = i\sin\theta$$

can be used in the case *c*<sup>0</sup> < 0 to obtain

$$P \approx P\_0 \left( \cos \left( \sqrt{|c\_0|} \,\Delta r \right) - \frac{d\_0}{\sqrt{|c\_0|}} \sin \left( \sqrt{|c\_0|} \,\Delta r \right) \right),$$

which is useful whenever there is a computational issue with imaginary numbers.

In the next section, we demonstrate that the hyperbolic approximation is less prone to errors than other well-known approximations in situations where the duration and/or convexity are negative. Recall that negative convexity cash flow streams can be easily constructed with the addition of negative cash flows to a stream of positive payments, when considering callable bonds, or with mortgage-backed securities due to the prepayment option in conventional residential mortgages.

#### **3. Results**

Approximation formulas such as Equations (1)–(5) should ideally be intuitive and behave well in special cases.

(i) The simplest cash flow is cash, which has trivial duration and convexity and is unaffected by interest rate changes. By substituting *d*<sup>0</sup> = *c*<sup>0</sup> = 0 or taking the corresponding limit in the case of (5) and using the fact that

$$\lim\_{\theta \to 0} \frac{\sinh \theta}{\theta} = \lim\_{\theta \to 0} \frac{\sin \theta}{\theta} = 1,$$

we obtain *P* = *P*<sup>0</sup> as expected.


We supplement the theoretical tests above with some concrete examples.

(iv) Consider a 10-year annuity-immediate with annual payments of 10. Recall that our assumption is *r*<sup>0</sup> = 1.6% and compute the present value *P*<sup>0</sup> = 10*a* <sup>10</sup> = 91.6728. Another easy calculation gives the Macaulay duration and convexity as *d*<sup>0</sup> = 5.3681 and *c*<sup>0</sup> = 37.0554, respectively.

In Table 1, the exact value of *P* is computed using the same formula as for *P*<sup>0</sup> but at the new continuously compounded rate.


**Table 1.** PV of 10-year annuity with annual payments of 10.

We can observe that Tchuindjo's approximation outperforms the rest, while Barber's lags behind for sizable rate changes. This was to be expected, since Barber's approximation lacks a convexity term and will not do well in cases when *c*<sup>0</sup> − *d* 2 0 is non-trivial. On the other hand, Alps' and the hyperbolic approximations are roughly equally accurate behind Tchuindjo's.

(v) Next, add a negative cash flow at time 20. We have chosen *CF*<sup>20</sup> = −120 in the example below; the net present value is *P*<sup>0</sup> = 4.5349 and the Macaulay duration and convexity are *d*<sup>0</sup> = −275.7817 and *c*<sup>0</sup> = −6, 936.8498, respectively.

Looking at Table 2 below, it may come as a surprise that the approximations by Tchuindjo and Alps blow up completely. However, we can provide a simple mathematical explanation for the bizarre behavior. Whenever *c*<sup>0</sup> < 0, the quadratic term of these two approximations that includes the expression *c*<sup>0</sup> − *d* 2 0 has the potential to be extremely influential. As ∆*r* increases, said term can overwhelm the baseline value *P*<sup>0</sup> and the linear term, resulting in large errors. Barber's approximation exhibits the opposite weakness: missing a quadratic term implies that the negative convexity is not accounted for at all. In fact, for suitable *CF*20, we can obtain *d*<sup>0</sup> = 0, in which case Barber's approximation fails to yield any results.


**Table 2.** NPV of 10-year annuity with annual payments of 10 and a payment of −120 at time 20.

We conclude this example by mentioning that the top-performing approximation is Fischer–Weil's, while the hyperbolic approximation is second-best.

(vi) Let us now consider a dividend stock, whose theoretical price is computed using *Gordon's dividend discount model*

$$P = \frac{D}{r - g}$$

with *D* representing next year's dividend and *g* its constant continuously compounded growth rate in perpetuity. A quick calculation gives *d* = (*r* − *g*) <sup>−</sup><sup>1</sup> and *<sup>c</sup>* <sup>=</sup> <sup>2</sup>(*<sup>r</sup>* <sup>−</sup> *<sup>g</sup>*) −2 ; for *g* = 0.6% we obtain *d*<sup>0</sup> = 100 and *c*<sup>0</sup> = 20, 000. Assume *D* = 1.

Some of the results in Table 3 may appear counterintuitive at first sight. Gordon's model suggests that *P* has an inverse relationship to *r*; however, all approximations except for Alps' and Barber's eventually produce a divergent estimate for *P* as ∆*r* increases. However, this is explained by the fact that we are attempting to trace a hyperbola using quadratic curves. Moreover, all approximations struggle to keep up with *P* for large negative values of ∆*r*.

**Table 3.** Price of a dividend stock with *D* = 1 and *g* = 0.6%.


Overall, Alps' approximation proves to be the most dependable for moderate changes in the interest rates.

(vii) Next, consider a 10-year bond with a coupon rate of *r*<sup>0</sup> and face value of 100. A quick calculation yields *d*<sup>0</sup> = 9.3151 and *c*<sup>0</sup> = 90.6932.

It turns out that the last three approximation formulas clearly outperform the rest, with Tchuindjo's having a slight advantage over Alps' and the hyperbolic approximation, as evidenced from Table 4. The subpar performance of Fischer–Weil on bonds is one of the reasons why this approximation is not widely utilized, despite its robustness in cases such as (v).

It has been shown empirically that although investment-grade bonds fall in price when interest rates rise, that is not necessarily the case with high-yield bonds whose duration can be negative due to default risk; see Melentyev and Yu (2020). For such bonds, care should be exercised when using the approximations by Tchuindjo or Alps.


**Table 4.** PV of 10-year par bond with coupon rate *r*0.

(viii) Finally, assume the bond is callable, with the European call strike set at *K* = 101.0000 and bond price volatility *σ* = 8%. The call is exercised a year ahead of the bond's maturity and has price *V* = 9.9431, which is subtracted from the price of a conventional bond to arrive at the callable bond price. Using ∆*r* = 20 bp in the calculation of the effective duration and convexity, we obtain *d<sup>e</sup>* = 5.3333 and *c<sup>e</sup>* = 50.3993. The positive convexity may surprise some readers, but note that the convexity turns negative when the interest rate gets closer to 0 and the bond price approaches the strike.

It is important to observe in Table 5 that none of the approximation formulas can consistently outperform the rest, if our objective is to estimate the full range of prices for such a bond. In effect, we are trying to approximate a function with an inflection point using quadratic curves, and thus significant approximation errors are inevitable.


**Table 5.** NPV of 10-year par bond with coupon rate *r*0, callable for 101 at *t* = 9.

The only useful conclusion is that the hyperbolic approximation is never the worst one, since it tends to be "sandwiched" between other approximations.

#### **4. Discussion**

We have established through a number of theoretical considerations and concrete examples that the accuracy of various Macaulay approximations can vary widely. Approximations that outperform in one case turn out to be unreliable in another case. The hyperbolic approximation, introduced in this paper, exhibited modest errors in most cases and thus the most reliability among the five approximations studied.

We can envision a variety of uses for the results presented here:

• To perform expeditious interest risk calculations by practitioners;


There is also potential to expand the scope of this study by incorporating non-flat yield curves, key rate durations, passage of time, and more complex financial instruments.

**Funding:** This research received no external funding.

**Acknowledgments:** The author wishes to thank the anonymous referees for their constructive comments. Partial results were presented at the 56th Actuarial Research Conference (see Orfanos 2021) and this work benefitted from the attendees' questions.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


Tchuindjo, Leonard. 2008. An Accurate Formula for Bond-Portfolio Stress Testing. *The Journal of Risk Finance* 9: 262–77. [CrossRef]

## *Article* **Reverse Sensitivity Analysis for Risk Modelling**

**Silvana M. Pesenti**

Department of Statistical Sciences, University of Toronto, Toronto, ON M5S 3G3, Canada; silvana.pesenti@utoronto.ca

**Abstract:** We consider the problem where a modeller conducts sensitivity analysis of a model consisting of random input factors, a corresponding random output of interest, and a baseline probability measure. The modeller seeks to understand how the model (the distribution of the input factors as well as the output) changes under a stress on the output's distribution. Specifically, for a stress on the output random variable, we derive the unique stressed distribution of the output that is closest in the Wasserstein distance to the baseline output's distribution and satisfies the stress. We further derive the stressed model, including the stressed distribution of the inputs, which can be calculated in a numerically efficient way from a set of baseline Monte Carlo samples and which is implemented in the R package SWIM on CRAN. The proposed reverse sensitivity analysis framework is model-free and allows for stresses on the output such as (a) the mean and variance, (b) any distortion risk measure including the Value-at-Risk and Expected-Shortfall, and (c) expected utility type constraints, thus making the reverse sensitivity analysis framework suitable for risk models.

**Keywords:** distortion risk measures; expected utility; Wasserstein distance; robustness and sensitivity analysis; model uncertainty

#### **1. Introduction**

Sensitivity analysis is indispensable for model building, model interpretation, and model validation, as it provides insight into the relationship between model inputs and outputs. A key tool used for sensitivity analysis are sensitivity measures, that assign to each model input a score, representing an input factor's ability to explain the variability of a model output's summary statistic; see Saltelli et al. (2008) and Borgonovo and Plischke (2016) for an in-depth review. One of the most widely used output summary statistic is the variance, which gives rise to sensitivity measures, e.g., the Sobol indices, that apportion the uncertainty in the output's variance to input factors. In many applications, such as reliability management and financial and insurance risk management, however, the variance is not the output statistic of concern and instead quantile-base measures are used; indicatively, see Asimit et al. (2019); Fissler and Pesenti (2022); Maume-Deschamps and Niang (2018); Tsanakas and Millossovich (2016). Furthermore, typical for financial risk management applications is that model inputs are subject to distributional uncertainty. Probabilistic (or global) sensitivity measures, however, tacitly assume that the model's distributional assumptions are correctly specified; indeed, sensitivity measures based on the difference between conditional (on a model input) and unconditional densities (of the output) are termed "common rationale" Borgonovo et al. (2016). Examples include indices, such as Borgonovo's sensitivity measures Borgonovo (2007), the *f*-sensitivity index Rahman (2016), and sensitivity indices based on the Cramér–von Mises distance Gamboa et al. (2018), we also refer to Plischke and Borgonovo (2019) for a detailed overview and to Gamboa et al. (2020) for estimation of these sensitivity measures. Recently, Plischke and Borgonovo (2019) define sensitivity measures that depend only on the copula between input factors, whereas Pesenti et al. (2021) propose a sensitivity measure based on directional derivatives that take dependence between input factors into account. Estimating these sensitivities, however, may render difficult in application where joint observations are

**Citation:** Pesenti, Silvana M. 2022. Reverse Sensitivity Analysis for Risk Modelling. *Risks* 10: 141. https:// doi.org/10.3390/risks10070141

Academic Editors: Annamaria Olivieri and Ermanno Pitacco

Received: 26 May 2022 Accepted: 7 July 2022 Published: 18 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

scarce, e.g., insurance portfolios, and their interpretation may be limited as dependence structures are commonly specified by expert opinions Denuit et al. (2006).

We consider an alternative sensitivity analysis framework proposed in Pesenti et al. (2019) that (a) considers statistical summaries relevant to risk management, (b) applies to models subject to distributional uncertainty, thus instead of relying on correctly specified distributions from which to calculate sensitivity measures we derive alternative distributions that fulfil a specific probabilistic stress and are "closest" to the baseline distribution; and (c) studies reverse sensitivity measures. Differently to the framework proposed in Pesenti et al. (2019) who use the Kullback–Leibler divergence to quantify the closedness of probability measures, in this work we consider the Wasserstein distance of order two to measure the distance between distribution functions. The Wasserstein distance allows for more flexibility in the choice of stresses including survival probabilities (via quantiles) used in reliability analysis, risk measures employed in finance and insurance, and utility functions relevant for decision under ambiguity.

Central to the reverse sensitivity analysis framework is a *baseline model*, the 3-tuple (*X*, *g*, P), consisting of random input factors *X* = (*X*1, . . . , *Xn*), an aggregation function *<sup>g</sup>* : <sup>R</sup>*<sup>n</sup>* <sup>→</sup> <sup>R</sup> mapping input factors to a univariate output *<sup>Y</sup>* <sup>=</sup> *<sup>g</sup>*(*X*), and a probability measure P. The methodology has been termed *reverse sensitivity analysis* by Pesenti et al. (2019) since it proceeds in a reverse fashion to classical sensitivity analysis where input factors are perturbed and the corresponding altered output is studied. Indeed, in the reverse sensitivity analysis proposed by Pesenti et al. (2019) a stress on the output's distribution is defined and changes in the input factors are monitored. The quintessence of the sensitivity analysis methodology is, however, not confined to stressing the output's distribution, it is also applicable to stressing an input factor and observing the changes in the model output and in the other inputs. Throughout the exposition, we focus on the reverse sensitivity analysis that proceeds via the following steps:


Sensitivity testing using divergence measures–in the spirit of the reverse sensitivity methodology–has been studied by Cambou and Filipovi´c (2017) using *f*-divergences on a finite probability space; by Pesenti et al. (2019) and Pesenti et al. (2021) using the Kullback– Leibler divergence; and Makam et al. (2021) consider a discrete sample space combined with the *χ* 2 -divergence. It is however known that the set of distribution functions with finite *f*-divergence, e.g., the Kullback–Leibler and *χ* <sup>2</sup> divergence–around a baseline distribution function depends on the baseline's tail-behaviour, thus the choice of *f*-divergence should be chosen dependent on the baseline distribution Kruse et al. (2019). The Wasserstein distance on the contrary, automatically adapts to the baseline distribution function in that the Wasserstein distance penalises dissimilar distributional features such as different tail behaviour Bernard et al. (2020). The Wasserstein distance has enjoyed numerous applications to quantify distributional uncertainty, see, e.g., Blanchet and Murthy (2019) and Bernard et al. (2020) for applications to financial risk management. In the context of uncertainty quantification, Moosmüeller et al. (2020) utilise the Wasserstein distance to elicit the (uncertain) aggregation map *g* from the distributional knowledge of the inputs and outputs. Fort et al. (2021) utilises the Wasserstein distance to introduce global sensitivity indices for computer codes whose output is a distribution function. In this manuscript we use the Wasserstein distance as it allows for different stresses compared to the Kullback–Leibler divergence. Indeed, the Wasserstein distance allows for stresses on any distortion risk measures, while the Kullback–Leibler divergence only allow for stresses on risk measures which are Value-at-Risk (VaR) and VaR and Expected Shortfall jointly, see Pesenti et al. (2019).

This paper is structured as follows: In Section 2, we state the notation and definitions necessary for the exposition. Section 3 introduces the optimisation problems and we derive the unique stressed distribution function of the output which has minimal Wasserstein distance to the baseline output's distribution and satisfies a stress. The considered stresses include constraints on risk measures, quantiles, expected utilities, and combinations thereof. In Section 4, we characterise the canonical Radon–Nikodym derivative, induced by the stressed distribution function, and study how input factors' distributions change when moving from the baseline to the stressed model. An application of the reverse sensitivity analysis is demonstrated on a mixture model in Section 5.

All proofs are delegated to Appendix A.

#### **2. Preliminaries**

Throughout we work on a measurable space (Ω, A) and denote the sets of distribution functions with finite second moment by

$$\mathcal{M} = \left\{ G \colon \mathbb{R} \to [0, 1] \; \middle| \; G \text{ non-decreasing} \text{ right-continuous} \; \begin{aligned} & \lim\_{\mathbf{x} \searrow -\infty} \mathcal{G}(\mathbf{x}) = 0, \\ & \lim\_{\mathbf{x} \nearrow +\infty} \mathcal{G}(\mathbf{x}) = 1, \text{ and } \int \mathbf{x}^2 d\mathcal{G}(\mathbf{x}) < +\infty \end{aligned} \right\},$$

and the corresponding set of square-integrable (left-continuous) quantile functions by

$$\mathcal{M} = \left\{ \begin{array}{c} \mathsf{\tilde{G}} \in \mathbb{L}^{2}([0,1]) \mid \mathsf{\tilde{G}} \text{ non-decreasing} \\ \end{array} \; \middle| \; \begin{array}{c} \mathsf{left\prime} \text{non-decreasing} \\ \end{array} \right\}.$$

or any distribution function *G* ∈ M, we denote its corresponding (left-continuous) quantile function by *<sup>G</sup>*˘ <sup>∈</sup> <sup>M</sup>˘ , that is *<sup>G</sup>*˘(*u*) = inf{ *<sup>y</sup>* <sup>∈</sup> <sup>R</sup> <sup>|</sup> *<sup>G</sup>*(*y*) <sup>≥</sup> *<sup>u</sup>*}, *<sup>u</sup>* <sup>∈</sup> [0, 1], with the convention that inf ∅ = +∞. We measure the discrepancy between distribution functions on the real line using the Wasserstein distance of order 2, defined as follows.

**Definition 1** (Wasserstein Distance)**.** *The Wasserstein distance (of order 2) between two distribution functions F*<sup>1</sup> *and F*<sup>2</sup> *is defined as Villani (2008)*

$$\mathcal{W}\_2(F\_1, F\_2) = \inf\_{\pi \in \Pi(F\_1, F\_2)} \left\{ \left( \int\_{\mathbb{R}^2} |z\_1 - z\_2|^2 \, \pi(dz\_1, dz\_2) \right)^{\frac{1}{2}} \right\},$$

*where* Π(*F*1, *F*2) *denotes the set of all bivariate probability measures with marginal distributions F*<sup>1</sup> *and F*2*, respectively.*

The Wasserstein distance is the minimal quadratic cost associated with transporting the distribution *F*<sup>1</sup> to *F*<sup>2</sup> using all possible couplings (bivariate distributions) with fixed marginals *F*<sup>1</sup> and *F*2. The Wasserstein distance admits desirable properties to quantify model uncertainty such as the comparison of distributions with differing support, e.g., with the empirical distribution function. Moreover it is symmetric and forms a metric on the space of probability measures; we refer to Villani (2008) for an overview and properties of the Wasserstein distance. It is well known (Dall'Aglio 1956) that for distributions on the real line, the Wasserstein distance admits the representation

$$\mathcal{W}\_2(F\_1 \mid F\_2) = \left( \int\_0^1 \left| \check{F}\_1(u) - \check{F}\_2(u) \right|^2 du \right)^{\frac{1}{2}}.$$

#### **3. Deriving the Stressed Distribution**

Throughout this section we assume that the modeller's *baseline model* is the 3-tuple (*X*, *g*, P) consisting of a random vector of input factors *X* = (*X*1, . . . , *Xn*), an aggregation function *<sup>g</sup>* : <sup>R</sup>*<sup>n</sup>* <sup>→</sup> <sup>R</sup> mapping input factors to a (for simplicity) univariate output *<sup>Y</sup>* <sup>=</sup> *<sup>g</sup>*(*X*), and a probability measure P. The baseline probability measure P reflects the modeller's (statistical and expert) knowledge of the distribution of *X* and we denote the distribution function of the output by *F*(*y*) = P(*Y* ≤ *y*). The modeller then performs reverse sensitivity analysis, that is tries to understand how prespecified stresses/constraints on the output distribution *F*, e.g., an increase in jointly its mean and standard deviation or a risk measures such as the Value-at-Risk (VaR) or Expected Shortfall (ES), affects the baseline model, e.g., the joint distribution of the input factors. For this, we first define the notion of a *stressed distribution*. Specifically, for given constraints we call a solution to the optimisation problem

$$\underset{G \in \mathcal{M}}{\text{arg min }} \mathcal{W}\_2(G, F) \qquad \text{subject to stresses/constraints on } G \text{ } \tag{1}$$

a stressed distribution. In problem (1), the baseline distribution *F* is fixed and we seek over all alternative distributions *G* ∈ M the one who satisfies the stress(es) and which has smallest Wasserstein distance to *F*. The solution to problem (1)–the stressed distribution– may be interpreted as the most "plausible" distribution function arising under adverse circumstances. Examples of stresses and constraints considered in this work include an increase (decrease), compared to their corresponding values under the reference probability P, in e.g., the mean, mean and standard deviation, distortion risk measures, and utility functions, and combinations thereof.

Next, we recall the concept of weighted isotonic projection which is intrinsically connected to the solution of optimisation problem (1); indeed the stressed quantile functions can be uniquely characterised via weighted isotonic projections.

**Definition 2** (Weighted Isotonic Projection Barlow et al. (1972))**.** *The weighted isotonic projection* ` <sup>↑</sup>*<sup>w</sup> of a function* ` ∈ L 2 ([0, 1]) *with weight function w*: [0, 1] → [0, +∞)*, w* ∈ L 2 ([0, 1])*, is its weighted projection onto the set of non-decreasing and left-continuous functions in* L 2 ([0, 1])*. That is, the unique function satisfying*

$$\ell^{\uparrow\_w} = \underset{\hbar \in \breve{\mathcal{M}}}{\text{arg min }} \int\_0^1 \left( \ell(\mu) - h(\mu) \right)^2 w(\mu) \, d\mu.$$

When the weight function is constant, i.e., *w*(*x*) ≡ *c*, *c* > 0, we write ` ↑ (·) = ` <sup>↑</sup>*<sup>c</sup>* (·), as in this case the isotonic projection is indeed independent of *c*. The weighted isotonic projection admits not only a graphical interpretation as the non-decreasing function that minimises the weighted L 2 -distance from ` but has also a discrete counterpart: the weighted isotonic regression Barlow et al. (1972). Numerically efficient algorithms for calculating weighted isotonic regressions are available, e.g., the R package isotone De Leeuw et al. (2010).

In the next sections, we solve problem (1) for different choices of constraints. Specifically, for risk measures constraints (Section 3.1), integral constraints (Section 3.2), Value-at-Risk constraints (Section 3.3), and expected utility constraint (Section 3.4), and in Section 3.5 we consider ways to smooth stressed distributions. Using these stressed distributions, we derive the stressed probability measures in Section 4 and study how a stress on the output is reflected on the input distribution(s).

#### *3.1. Risk Measure Constraints*

This section considers stresses on distortion risk measures, that is we derive the unique stressed distribution that satisfies an increase and/or decrease of distortion risk measures while minimising the Wasserstein distance to the baseline distribution *F*.

**Definition 3** (Distortion Risk Measures)**.** *Let γ* ∈ L 2 ([0, 1]) *be a square-integrable function with <sup>γ</sup>*: [0, 1] <sup>→</sup> [0, <sup>+</sup>∞) *and* <sup>R</sup> <sup>1</sup> 0 *γ*(*u*) *du* = 1*. Then the distortion risk measure ρ<sup>γ</sup> with distortion weight function γ is defined as*

$$\rho\_{\mathcal{I}}(\mathcal{G}) = \int\_0^1 \mathcal{G}(u)\gamma(u) \, du \quad \text{for} \quad \mathcal{G} \in \mathcal{M}. \tag{2}$$

The above definition of distortion risk measures makes the assumption that positive realisations are undesirable (losses) while negative realisations are desirable (gains). The class of distortion risk measures includes one of the most widely used risk measures in financial risk management, the Expected Shortfall (ES) at level *α* ∈ [0, 1) (also called Tail Value-at-Risk), with *γ*(*u*) = <sup>1</sup> 1−*α* <sup>1</sup>{*u*>*α*} , see, e.g., Acerbi and Tasche (2002). The often used risk measure Value-at-Risk (VaR), while admitting a representation given in (2), has a corresponding weight function *γ* that is not square-integrable. We derive the solution to optimisation problem (1) with a VaR constraint in Section 3.3.

**Theorem 1** (Distortion Risk Measures)**.** *Let r<sup>k</sup>* ∈ R*, ργ<sup>k</sup> be a distortion risk measure with weight function <sup>γ</sup><sup>k</sup> and assume there exists a distribution function <sup>G</sup>*˜ ∈ M *satisfying <sup>ρ</sup>γ<sup>k</sup>* (*G*˜) = *r<sup>k</sup> for all k* ∈ {1, . . . , *d*}*. Then, the optimisation problem*

$$\underset{\mathbf{G}\in\mathcal{M}}{\arg\min} \ W\_2(\mathbf{G}, \mathbf{F}) \qquad \text{subject to} \qquad \rho\_{\gamma\_k}(\mathbf{G}) = r\_k \quad k = 1, \dots, d,\tag{3}$$

*has a unique solution given by*

$$\mathbf{G}^\*(\boldsymbol{\mu}) = \left(\check{\mathbf{F}}(\boldsymbol{\mu}) + \sum\_{k=1}^d \lambda\_k \gamma\_k(\boldsymbol{\mu})\right)^\uparrow,\tag{4}$$

*where the Lagrange multipliers λ<sup>k</sup> are such that the constraints are fulfilled, that is ργ<sup>k</sup>* (*G* ∗ ) = *r<sup>k</sup> for all k* = 1, . . . , *d.*

In the above theorem, and also in later results, we assume that there exists a distribution function which satisfies all constraints. This assumption is not restrictive and requires that, particularly, multiple constraints are chosen carefully, e.g., imposing that R 1 0 *G*˘(*u*)*du* > <sup>1</sup> 1−*α* R 1 *α <sup>G</sup>*˘(*u*)*du* for *<sup>α</sup>* <sup>∈</sup> (0, 1), i.e., the mean being larger than the ES*α*, cannot be fulfilled by any distribution function; thus, a combination of stresses not of interest to a modeller.

We observe that the optimal quantile function is the isotonic projection of a weighted linear combination of the baseline's quantile function *F*˘ and the distortion weight functions of the risk measures. A prominent group of risk measures is the class of coherent risk measures, that are risk measures fulfilling the properties of monotonicity, positive homogeneity, translation invariance, and sub-additivity; see Artzner et al. (1999) for a discussion and interpretation. It is well-known that a distortion risk measure is coherent, if and only if, its distortion weight function *γ*(·) is non-decreasing Kusuoka (2001). For the special case of a constraint on a coherent distortion risk measure that results in a larger risk measure compared to the baseline's, we obtain an analytical solution without the need to calculate an isotonic projection.

**Proposition 1** (Coherent Distortion Risk Measure)**.** *If ρ<sup>γ</sup> is a coherent distortion risk measure and r* ≥ *ργ*(*F*)*, then optimisation problem* (3) *with d* = 1 *has a unique solution given by*

$$\check{G}^\*(u) = \check{F}(u) + \frac{r - \rho\_\gamma(F)}{\int\_0^1 (\gamma(u))^2 du} \,\gamma(u) \,\,\iota$$

We illustrate the stressed distribution functions for constraints on distortion risk measures in the next example. Specifically, we look at the *α*-*β* risk measures which are a parametric family of distortion risk measures.

**Example 1** (*α*-*β* Risk Measure)**.** *The α-β risk measure,* 0 < *β* ≤ *α* < 1*, is defined by*

$$\gamma(\mu) = \frac{1}{\eta} \left( p \operatorname{\mathbf{1}}\_{\{\mu < \beta\}} + (1 - p) \operatorname{\mathbf{1}}\_{\{\mu \ge \alpha\}} \right),$$

*where p* ∈ [0, 1] *and η* = *p β* + (1 − *p*) (1 − *α*) *is the normalising constant. This parametric family contains several notable risk measures as special cases: for p* = 0 *we obtain ESα, and for p* = 1 *the conditional lower tail expectation (LTE) at level β.*

*Moreover, if p* < <sup>1</sup> 2 *p* > <sup>1</sup> 2 *the α-β risk measure emphasises losses (gains) relative to gains (losses). For α* = *β and p* < <sup>1</sup> 2 *, the risk measure is equivalent to κ*(*ESα*[*Y*] − *λ* E[*Y*])*, where κ* = (1−2*p*) (1−*α*) *η and λ* = *p κ η .*

*Figure 1 displays the baseline F*˘ *<sup>Y</sup> and the stressed G*˘ <sup>∗</sup> *Y quantile functions of a random variable Y under a 10% increase on the α-β risk measure with β* = 0.1*, α* = 0.9*, and various p* ∈ {0.25, 0.5, 0.75}*. The baseline distribution is chosen to be F<sup>Y</sup> is Lognormal*(*µ*, *σ* 2 ) *with parameters µ* = <sup>7</sup> 8 *and σ* = 0.5*. We observe in Figure 1 that the stressed quantile functions G*˘ <sup>∗</sup> *Y have, in all three plots, a flat part which straddles β* = 0.1 *and a jump at α* = 0.9*. The length of the flat part is increasing with increasing p while the size of the jump is decreasing with increasing p. This can also be seen in the stressed densities g* ∗ *<sup>Y</sup> which have, for all values of p, a much heavier right albeit a much lighter left tail than the density of the baseline model. Thus, under this stress, both tails of the baseline distribution are altered.*

**Figure 1. Top** panels: Baseline quantile function *F*˘ *<sup>Y</sup>* (blue dashed) compared to the stressed quantile function *G*˘ <sup>∗</sup> *Y* (red solid) for a 10% increase on the *α*-*β* risk measure with *β* = 0.1, *α* = 0.9, and various values of *<sup>p</sup>*. The green line `(·) is the function, whose isotonic projection equals *<sup>G</sup>*˘ *<sup>Y</sup>*(·). **Bottom** panels: corresponding baseline *f<sup>Y</sup>* and stressed *g* ∗ *Y* densities.

#### *3.2. Integral Constraints*

The next results are generalisations of stresses on distortion risk measures to integral constraints, and include as a special case a stress jointly on the mean, the variance, and distortion risk measures.

**Theorem 2** (Integral)**.** *Let h<sup>k</sup>* , ˜*hl* : [0, 1] → [0, ∞) *be square-integrable functions and assume there exists a distribution function <sup>G</sup>*˜ ∈ M *satisfying* <sup>R</sup> <sup>1</sup> 0 *hk* (*u*)*G*˘(*u*) *du* <sup>≤</sup> *<sup>c</sup><sup>k</sup> and* R 1 0 ˜*hl* (*u*) *G*˘(*u*) 2 *du* ≤ *c*˜*<sup>l</sup> for all k* = 1, . . . , *d, and l* = 1, . . . , ˜*d. Then the optimisation problem*

$$\begin{aligned} \underset{\mathbf{G}\in\mathcal{M}}{\operatorname{arg\,min}}\ W\_2(\mathbf{G},\mathbf{F}) \qquad \text{subject to} \qquad \int\_0^1 h\_k(u)\check{\mathbf{G}}(u)\,du \le \mathfrak{c}\_k, \quad k=1,\dots,d, \\\int\_0^1 \tilde{h}\_l(u) \left(\check{\mathbf{G}}(u)\right)^2 du \le \mathfrak{c}\_l, \quad l=1,\dots,\tilde{d}\_l. \end{aligned}$$

*has a unique solution given by*

$$\mathcal{G}^\*(u) = \left(\frac{1}{\tilde{\Lambda}(u)} \left(\tilde{\mathcal{F}}(u) + \sum\_{k=1}^d \lambda\_k h\_k(u)\right)\right)^{\uparrow\_{\tilde{\Lambda}}},$$

*where* Λ˜ (*u*) = 1 + ∑ ˜*d k*=1 *λ*˜ *k* ˜*hk* (*u*) *and the Lagrange multipliers λ*1, . . . , *λ<sup>d</sup> and λ*˜ <sup>1</sup>, . . . , *λ*˜ *<sup>d</sup> are non-negative and such that the constraints are fulfilled.*

A combination of the above theorems provides stresses jointly on the mean, the variance, and on multiple distortion risk measures.

**Proposition 2** (Mean, Variance, and Risk Measures)**.** *Let m*<sup>0</sup> ∈ R*, σ* <sup>0</sup> > 0*, r<sup>k</sup>* ∈ R*, and distortion risk measures ργ<sup>k</sup> , <sup>k</sup>* <sup>=</sup> 1, . . . , *d. Assume there exists a distribution function <sup>G</sup>*˜ ∈ M *with mean m*0 *, standard deviation σ* 0 *, and which satisfies ργ<sup>k</sup>* (*G*˜) = *r<sup>k</sup> , for all k* = 1, . . . , *d. Then the optimisation problem*

$$\begin{aligned} \underset{G \in \mathcal{M}}{\text{arg min }} \; W\_2(G, F) \qquad & \text{subject to} \qquad \int \mathbf{x} \, dG(\mathbf{x}) = m', \\ & \int (\mathbf{x} - \mathbf{m'})^2 \, dG(\mathbf{x}) = \left(\sigma'\right)^2 \quad \text{and} \\ & \rho\_{\mathcal{M}}(\mathbf{G}) = r\_{k'} \quad k = 1, \dots, d\_{\prime} \end{aligned}$$

*has a unique solution given by*

$$\mathbf{G}^\*(u) = \left(\frac{1}{1+\lambda\_2} \left(\check{\mathbf{F}}(u) + \lambda\_1 + \lambda\_2 m' + \sum\_{k=1}^d \lambda\_{k+2} \gamma\_k(u)\right)\right)^\uparrow \mathbf{f}$$

*and the Lagrange multipliers λ*1, . . . , *λd*+<sup>2</sup> *with λ*<sup>2</sup> 6= −1 *are such that the constraints are fulfilled.*

**Example 2** (Mean, Variance, and ES)**.** *Here, we illustrate Proposition 2 with the ES risk measure and three different stresses. The top panels of Figure 2 display the baseline quantile function F*˘ *<sup>Y</sup> and the stressed quantile function G*˘ <sup>∗</sup> *Y of Y, where the baseline distribution F<sup>Y</sup> of Y is again Lognormal*(*µ*, *σ* 2 ) *with parameters µ* = <sup>7</sup> 8 *and σ* = 0.5*. The bottom panels display the corresponding baseline and stressed densities. The left panels correspond to a stress, where, under the stressed model, the ES*0.95 *and the mean are kept fixed at their corresponding values under the baseline model, while the standard deviation is increased by 20%. We observe, both in the quantile and density plot, that the stressed distribution is more spread out indicating a larger variance. Furthermore, at y* ≈ 5.77 *the stressed density g* ∗ *Y* (*y*) *drops to ensure that ES*0.95(*G* ∗ *Y* ) = *ES*0.95(*FY*)*. This drop is due to the fact that a stress composed of a 20% increase in the standard deviation while fixing the mean (i.e., without a constraint on ES) results in an ES that is larger compared to the baseline's. Indeed, under this alternative stress (without a constraint on ES) we obtain that ES*0.95(*G* ∗ *Y* ) ≈ 7.70 *compared to ES*0.95(*FY*) ≈ 6.87*.*

*The middle panels correspond to a 10% increase in ES*0.95 *and a 10% decrease in the mean, while keeping the standard deviation fixed at its value under the baseline model. The density plot clearly indicates a general shift of the stressed density to the left, stemming from the decrease in the*

*mean, and a single trough which is induced by the increase in ES. The right panels correspond to a 10% increase in ES*0.95*, a 10% increase in the mean, and a 10% decrease in the standard deviation. The stressed density still has the trough from the increase in ES; however, the density is less spread out (reduction in the standard deviation) and generally shifted to the right (increase in the mean).*

**Figure 2. Top**: Baseline quantile function *F*˘ *<sup>Y</sup>* compared to the stressed quantile function *G*˘ <sup>∗</sup> *Y* . **Bottom**: corresponding baseline *f<sup>Y</sup>* and stressed *g* ∗ *Y* densities. **Left**: ES0.95 and the mean being fixed and a 20% increase in the standard deviation. **Middle**: 10% increase in ES0.95, 10% decrease in the mean, and fixed standard deviation. **Right**: 10% increase in ES0.95, 10% increase in the mean, and 10% decrease in standard deviation. Note that in the middle and right panel the green lines are equal to the red lines.

#### *3.3. Value-at-Risk Constraints*

In this section we study stresses on the risk measure Value-at-Risk (VaR). The VaR at level *α* ∈ (0, 1) of a distribution function *G* ∈ M is defined as its left-continuous quantile function evaluated at *α*, that is

$$\text{VaR}\_{\mathfrak{a}}(G) = \vec{G}(\mathfrak{a})\,.$$

We further define the right-continuous VaR+, that is the right-continuous quantile function of *G* ∈ M evaluated at *α*, by

$$\text{VaR}^+\_{\alpha}(G) = G^+(\alpha) = \inf \left\{ y \in \mathbb{R} \mid F(y) > \alpha \right\} \dots$$

**Theorem 3** (VaR)**.** *Let q* ∈ R *and consider the optimisation problem*

$$\begin{array}{ccc}\underset{G\in\mathcal{M}}{\arg\min} & W\_2(G, F) & \text{subject to} & (a) & VaR\_a(G) = q & or\\ & & (b) & VaR\_a^+(G) = q \end{array}$$

*and define α<sup>F</sup> such that VaRα<sup>F</sup>* (*F*) = *q. Then, the following holds*

*(i) under constraint (a), if q* ≤ *VaRα*(*F*)*, then the unique solution is given by*

$$\mathcal{G}^\*(\mu) = \mathring{F}(\mu) + \left(q - \mathring{F}(\mu)\right) \mathbb{1}\_{\{\mu \in (\mathfrak{a}\_F, \mathfrak{a}]\}} \grave{ }^\circ \check{\mathsf{F}}$$

*α*

*if q* > *VaRα*(*F*)*, then there does not exist a solution.*

*(ii) under constraint (b), if q* <sup>≥</sup> *VaR*<sup>+</sup> *α* (*F*)*, then the unique solution is given by*

$$\check{G}^\*(\mu) = \mathcal{F}(\mu) + \left(q - \mathcal{F}(\mu)\right) \mathbb{1}\_{\{\mu \in (\mathfrak{a}, \mathfrak{a}\_F]\}} \circ$$

*if q* < *VaR*<sup>+</sup> *α* (*F*)*, then there does not exist a solution.*

The above theorem states that if the optimal quantile function exists it is either the baseline quantile function *F*˘ or constant equal to *q*. Moreover, the stressed quantile function (if it exists) jumps at *α* which implies that the existence of a solution hinges on the careful choice of the stress. For a stress on VaR (constraint (a)) for example, a solution exists if and only if the constraint satisfies *q* ≤ VaR*α*(*F*); a decrease in the VaR*<sup>α</sup>* from the baseline to the stressed model. The reason for the non-existence of a solution when stressing VaR upwards is that the unique increasing function that minimises the Wasserstein distance and satisfies the constraint is not left-continuous and thus not a quantile function.

Alternatively to stressing VaR or VaR+, and in particularly in the case when a desired stressed solution does not exist, one may stress instead the distortion risk measure Range-Value-at-Risk (RVaR) Cont et al. (2010). The RVaR at levels 0 ≤ *α* < *β* ≤ 1 is defined by

$$\text{RVaR}\_{\mathfrak{n},\beta}(G) = \frac{1}{\beta - \mathfrak{a}} \int\_{\mathfrak{a}}^{\beta} \breve{G}(\mathfrak{n}) \, d\mathfrak{n} \, , \quad \text{for} \quad G \in \mathcal{M}\_{\omega}$$

and belongs to the class of distortion risk measures. The RVaR attains as limiting cases the VaR and VaR+. Indeed, for any *<sup>G</sup>* ∈ M it holds

$$\text{VaR}\_{\mathfrak{a}}(\mathcal{G}) = \lim\_{\mathfrak{a}' \nearrow \mathfrak{a}} \text{RVR}\_{\mathfrak{a}', \mathfrak{a}}(\mathcal{G}) \quad \text{and} \quad \text{VaR}\_{\mathfrak{a}}^+(\mathcal{G}) = \lim\_{\mathcal{F} \searrow \mathfrak{a}} \text{RVR}\_{\mathfrak{a}, \mathfrak{f}}(\mathcal{G}) \,.$$

The solution to stressing RVaR is provided in Theorem 1.

#### *3.4. Expected Utility Constraint*

This section considers the change from the baseline to the stressed distribution under an increase of an expected utility constraint. In the context of utility maximisation, the next theorem provides a way to construct stressed models with a larger utility compared to the baseline.

**Theorem 4** (Expected Utility and Risk Measures)**.** *Let u*: R → R *be a differentiable concave utility function, r<sup>k</sup>* ∈ R*, and ργ<sup>k</sup> be distortion risk measures, for k* = 1, . . . , *d. Assume there exists a distribution function G*˜ *satisfying* R R *<sup>u</sup>*(*x*) *dG*˜(*x*) <sup>≥</sup> *<sup>c</sup> and <sup>ρ</sup>γ<sup>k</sup> G*˜ = *r<sup>k</sup> for all k* = 1, . . . , *d. Then the optimisation problem*

$$\underset{\mathbf{G}\in\mathcal{M}}{\operatorname\*{arg\,min}}\ \,\,\mathcal{W}\_{2}(\mathbf{G},\mathbf{F})\qquad\text{subject to}\qquad\int\_{\mathbb{R}}\mathfrak{u}(\mathbf{x})\,\,d\mathbf{G}(\mathbf{x})\geq\mathbf{c}\qquad\&\quad\rho\_{\gamma\_{k}}(\mathbf{G})=r\_{k},\quad k=1,\dots,d\\\mathbb{R}\text{ is a vector of }\mathbf{g}\in\mathbb{R}^{d}\text{ if and only if }\quad\int\_{\mathbb{R}}\mathfrak{u}(\mathbf{x})\,\,d\mathbf{g}=r\_{k}\leq\mathbf{c}.\qquad\text{by the following conditions:}\\\int\_{\mathbb{R}}\mathfrak{u}(\mathbf{x})\,\,d\mathbf{g}=r\_{k}\leq\mathbf{c}.\qquad\text{by the following conditions:}\\\int\_{\mathbb{R}}\mathfrak{u}(\mathbf{x})\,\,d\mathbf{g}=r\_{k}\leq\mathbf{c}.\qquad\text{by the following conditions:}\\\int\_{\mathbb{R}}\mathfrak{u}(\mathbf{x})\,\,d\mathbf{g}=r\_{k}\leq\mathbf{c}.\qquad\text{by the}$$

*has a unique solution given by*

$$\mathcal{G}^\*(u) = \tilde{\nu}\_{\lambda\_1} \left( \left( \tilde{\mathbf{f}}(u) + \sum\_{k=1}^d \lambda\_{k+1} \gamma\_k(v) \right)^\gamma \right),\tag{5}$$

*where ν*˘*λ*<sup>1</sup> *is the left-inverse of νλ*<sup>1</sup> (*x*) = *x* − *λ*<sup>1</sup> *u* 0 (*x*)*, and λ*<sup>1</sup> ≥ 0*,* (*λ*2, . . . , *λd*+<sup>1</sup> ) <sup>∈</sup> <sup>R</sup>*<sup>d</sup> are such that the constraints are fulfilled.*

The utility function in Theorem 4 need not be monotone, indeed the theorem applies to any differentiable concave function, without the need of an utility interpretation. Moreover, Theorem 4 also applies to differentiable convex (disutility) functions *u*˜ and constraint R R *u*˜(*x*) *dG*(*x*) ≤ *c*; a situation of interest in insurance premium calculations. In this case, the solution is given by (5) with *u*(*x*) = −*u*˜(*x*).

**Example 3** (HARA Utility & ES)**.** *The Hyperbolic absolute risk aversion (HARA) utility function is defined by*

$$
\mu(\mathbf{x}) = \frac{1-\eta}{\eta} \left( \frac{a\mathbf{x}}{1-\eta} + b \right)^{\eta} \lambda
$$

*with parameters a* > 0*, ax* <sup>1</sup>−*<sup>η</sup>* <sup>+</sup> *<sup>b</sup>* <sup>&</sup>gt; <sup>0</sup>*, and where <sup>η</sup>* <sup>≤</sup> <sup>1</sup> *guarantees concavity.*

*We again choose the baseline distribution F<sup>Y</sup> of Y to be Lognormal*(*µ*, *σ* 2 ) *with µ* = <sup>7</sup> 8 *and σ* = 0.5 *and consider utility parameters a* = 1*, b* = 5*, and η* = 0.5*. Figure 3 displays the baseline and the stressed quantile functions F*˘ *<sup>Y</sup> and G*˘ <sup>∗</sup> *Y , respectively, for a combined stress on the HARA utility and on ES at levels 0.8 and 0.95. Specifically, for all three stresses we decreasing ES*0.8 *by 10% and increasing ES*0.95 *by 10% compared to their values under the baseline model. Moreover, the HARA utility is increased by 0%, 1%, and 3%, respectively, corresponding to the panels from the left to the right. The flat part in the stressed quantile function G*˘ <sup>∗</sup> (*u*) *around u* = 0.8*, visible in all top panels of Figure 3, is induced by the decrease in ES*0.8 *while the jump at u* = 0.95 *is due to the increase in ES*0.95*. From the left to the right panel in Figure 3, we observe that the larger the stress on the HARA utility, the more the stressed quantile function shifts away from the baseline quantile function F*˘ *Y.*

**Figure 3. Top** panels: Baseline quantile function *F*˘ *<sup>Y</sup>* compared to the stressed quantile function *G*˘ <sup>∗</sup> *Y* , for a 10% decrease in ES0.8, and a 10% increase in ES0.95, and, from left to right, a 0%, 1%, and 3% increase in the HARA utility, respectively. The function `(·) (solid green) is the function whose isotonic projection equals *G*˘ <sup>∗</sup> (·). **Bottom** panels: Corresponding baseline *f<sup>Y</sup>* and stressed *g* ∗ *Y* densities.

#### *3.5. Smoothing of the Stressed Distribution*

We observe that the stressed quantile functions derived in Section 3 generally contain jumps and/or flat parts even if the baseline distribution is absolutely continuous. In situation where this is not desirable, one may consider a smoothed version of the stressed distributions. For this, we recall that the isotonic regression, the discrete counterpart of the weighted isotonic projection, of a function ` evaluated at *u*1, . . . , *u<sup>n</sup>* with positive weights *w*1, . . . , *wn*, is the solution to

$$\min\_{\mu\_1, \dots, \mu\_n} \sum\_{i=1}^n (u\_i - \ell(u\_i))^2 w\_i \quad \text{subject to} \quad u\_i \le u\_j, \quad i \le j. \tag{6}$$

There are numerous efficient algorithms that solve (6) most notably the pool-adjacentviolators (PAV) algorithm Barlow et al. (1972). It is well-known that the solution to the

isotonic regression contains flat parts and jumps. A smoothed isotonic regression algorithm, termed smooth pool-adjacent-violators (SPAV) algorithm, using an L 2 regularisation was recently proposed by Sysoev and Burdakov (2019). Specifically, they consider

$$\min\_{\mu\_1, \dots, \mu\_n} \sum\_{i=1}^n (u\_i - \ell(u\_i))^2 w\_i + \sum\_{i=1}^n \zeta\_i (u\_{i+1} - u\_i)^2 \quad \text{subject to} \quad u\_i \le u\_j, \quad i \le j.$$

where *ζ<sup>i</sup>* ≥ 0, *i* = 0, . . . , *n* − 1, are prespecified smoothing parameters. Using a probabilitistic reasoning, Sysoev and Burdakov (2019) argue that *ζ<sup>i</sup>* may be chosen proportional to a (e.g., quadratic) kernel evaluated at *u<sup>i</sup>* and *ui*+<sup>1</sup> , that is

$$\mathcal{Z}\_i = \mathcal{Z}\,\mathsf{K}(\boldsymbol{u}\_i, \boldsymbol{u}\_{i+1}) \quad \text{with} \quad \mathsf{K}(\boldsymbol{u}\_i, \boldsymbol{u}\_{i+1}) = \frac{1}{|\boldsymbol{u}\_i - \boldsymbol{u}\_{i+1}|^2} \quad \text{and} \quad \mathcal{Z} \ge \mathbf{0} \,\mathsf{A}$$

The choice of smoothing parameter *ζ* = 0 correspond to the original isotonic regression larger values of *ζ* correspond to a greater degree of smoothness of the solution. *ζ* can either be prespecified or estimated using cross-validation, see e.g., Sysoev and Burdakov (2019).

To guarantee that the smoothed quantile function still fulfils the constraint, one may replace in every step of the optimisation for finding the Lagrange parameter the PAV with the SPAV algorithm. Thus, the Lagrange parameter are indeed found such that the constraints are fulfilled.

**Remark 1.** *There are numerous works proposing smooth versions of isotonic regressions. Approaches include kernel smoothers, e.g., Hall and Huang (2001), and spline techniques, e.g., Meyer (2008). These algorithms, however, are computationally heavy in that their computational cost is O*(*n* 2 )*, where n is the number of data points. Furthermore, these algorithm require a careful choice of the kernel or the spline basis which is in contrast to the SPAV. We refer the reader to Sysoev and Burdakov (2019) for a detailed discussion and references to smooth isotonic regression algorithms.*

#### **4. Analysing the Stressed Model**

Recall that a modeller is equipped with a baseline model, the 3-tuple (*X*, *g*, P), consisting of a set of input factors *X* = (*X*1, . . . , *Xn*), a univariate output random variable of interest, *Y* = *g*(*X*), and a probability measure P. For a stress on the output's baseline distribution *FY*, we derived in Section 3 the corresponding unique stressed distribution function, denoted here by *G* ∗ *Y* . Thus, to fully specify the stressed model we next define a stressed probability measure Q∗ that is induced by *G* ∗ *Y* .

#### *4.1. The Stressed Probability Measures*

A stressed distribution *G* ∗ *Y* induces a canonical change of measure that allows the modeller to understand how the baseline model including the distributions of the inputs changes under the stress. The Radon–Nikodym (RN) derivative of the baseline to the stressed model is

$$\frac{d\mathbb{Q}^\*}{d\mathbb{P}} = \frac{\mathbb{g}\_Y^\*(\mathcal{Y})}{f\_Y(\mathcal{Y})}$$

,

where *f<sup>Y</sup>* and *g* ∗ *Y* denote the densities of the baseline and stressed output distribution, respectively. The RN derivative is well-defined since *fY*(*Y*) > 0, P-a.s. The distribution functions of input factors under the stressed model – the stressed distributions – are then given, e.g., for input *X<sup>i</sup>* , *i* ∈ {1, . . . , *n*}, by

$$\mathbb{Q}^\*(X\_i \le \mathbf{x}\_i) = \mathbb{E}\left[\mathbf{1}\_{\{X\_i \le \mathbf{x}\_i\}} \frac{d\mathbb{Q}^\*}{d\mathbb{P}}\right] \\ = \mathbb{E}\left[\mathbf{1}\_{\{X\_i \le \mathbf{x}\_i\}} \frac{g\_Y^\*(Y)}{f\_Y(Y)}\right], \quad \mathbf{x}\_i \in \mathbb{R}\_{\ge 0}$$

and for multivariate inputs *X* by

$$\mathbb{Q}^\*(\mathbf{X} \le \mathbf{x}) = \mathbb{E}\left[\mathbf{1}\_{\{\mathbf{X} \le \mathbf{x}\}} \frac{\mathbf{g}\_Y^\*(Y)}{f\_Y(Y)}\right], \quad \mathbf{x} \in \mathbb{R}^n.$$

where E[·] denotes the expectation under P. Note that under the stressed probability measure Q∗ , the input factors' marginal and joint distributions may be altered.

**Example 4** (HARA Utility & ES continued)**.** *We continue Example 3 and illustrate the RNdensities <sup>d</sup>*Q<sup>∗</sup> *d*P *for the following three stresses (from the left to the right panel): a 10% decrease in ES*0.8 *and a 10% increasing ES*0.95 *for all three stresses, and a 0%, 1%, and 3% increase in the HARA utility, respectively.*

*We observe in Figure 4, that for all three stresses large realisations of Y obtain a larger weight under the stressed probability measures* Q∗ *compared to the baseline probability* P*. Indeed, for all three stresses it holds that <sup>d</sup>*Q<sup>∗</sup> *d*P (*ω*) > 1 *whenever Y*(*ω*) > 6 *and ω* ∈ Ω*. This is in contrast to small realisations of Y which obtain a weight smaller than 1. The impact of the different levels of stresses of the HARA utility (0%, 1%, and 3%, from the left to the right panel) can be observed in the left tail of <sup>d</sup>*Q<sup>∗</sup> *d*P *; a larger stress on the utility induces larger weights. The length of the trough of d*Q∗ *d*P *is increasing from the left panel (approx. (*4.53, 6.15*)) to the right panel (approx. (*4.43, 6.18*)), and corresponds in all cases to the constant part in G* ∗ *Y (see Figure 3, top panels) which is induced by the decrease in ES*0.8 *under the stressed model.*

**Figure 4.** RN-densities for following stresses: a 10% decrease in both ES0.8 and ES0.95, and an increase in the HARA utility. The change in HARA utility is 0%, 1%, and 3%, respectively, from left to right.

#### *4.2. Reverse Sensitivity Measures*

Comparison of the baseline and a stressed model can be conducted via different approaches depending on the modeller's interest. While probabilistic sensitivity measures underlie the assumption of a fixed probability measure and quantify the divergence between the conditional (on a model input) and the unconditional output density Saltelli et al. (2008), the proposed framework compares a baseline and a stressed model, i.e., distributions under different probability measures. Therefore, to quantify the distributional change in input factor *X<sup>i</sup>* from the baseline P to the stressed Q∗ probability, a sensitivity measure introduced by Pesenti et al. (2019) may be suitable which quantifies the variability of an input factor's distribution from the baseline to the stressed model. A generalisation of the reverse sensitivity measure is stated here.

**Definition 4** (Marginal Reverse Sensitivity Measure Pesenti et al. (2019))**.** *For a function s*: R → R*, the reverse sensitivity measure to input X<sup>i</sup> with respect to a stressed probability measure* Q∗ *is defined by*

$$\mathcal{S}\_{i}^{\mathcal{Q}^\*} = \begin{cases} \frac{\mathbb{E}^{\mathcal{Q}^\*}[s(X\_i)] - \mathbb{E}[s(X\_i)]}{\max\limits\_{\mathcal{Q} \in \mathcal{Q}} \mathbb{E}^{\mathcal{Q}}[s(X\_i)] - \mathbb{E}[s(X\_i)]} & \mathbb{E}^{\mathcal{Q}^\*}[s(X\_i)] \ge \mathbb{E}[s(X\_i)],\\ -\frac{\mathbb{E}^{\mathcal{Q}^\*}[s(X\_i)] - \mathbb{E}[s(X\_i)]}{\min\limits\_{\mathcal{Q} \in \mathcal{Q}} \mathbb{E}^{\mathcal{Q}}[s(X\_i)] - \mathbb{E}[s(X\_i)]} & \text{otherwise}, \end{cases}$$

*where* <sup>Q</sup> <sup>=</sup> {<sup>Q</sup> <sup>|</sup> <sup>Q</sup> *probability measure with <sup>d</sup>*<sup>Q</sup> *d*P P= *d*Q∗ *d*P } *is the set of all probability measures whose RN-derivative have the same distribution as <sup>d</sup>*Q<sup>∗</sup> *d*P *under* P*. We adopted the convention that* <sup>±</sup> <sup>∞</sup> <sup>∞</sup> <sup>=</sup> <sup>±</sup><sup>1</sup> *and* <sup>0</sup> <sup>0</sup> = 0*.*

The sensitivity measure is called "reverse", as the stress is applied to the output random variable *Y* and the sensitivity monitors the change in input *X<sup>i</sup>* . The definition of 4 applies, however, also to stresses on input factors, in which case the RN-density *d*Q∗ *d*P is a function of the stressed input factor and we refer to Pesenti et al. (2019) for a discussion. Note, that the reverse sensitivity measure can be viewed as a normalised covariance measure between the input *s*(*Xi*) and the Radon Nikodym derivative *<sup>d</sup>*Q<sup>∗</sup> *d*P .

The next proposition provides a collection of properties that the reverse sensitivity measure possesses, we also refer to Pesenti et al. (2019) for a detailed discussion of these properties. For this, we first recall the definition of comonotonic and counter-monotonic random variables.

**Definition 5.** *Two random variables Y*<sup>1</sup> *and Y*<sup>2</sup> *are comonotonic under* P*, if and only if, there exists a random variable W and non-decreasing functions h*1, *h*<sup>2</sup> : R → R*, such that the following equalities hold in distribution under* P

$$Y\_1 = h\_1(\mathcal{W}) \quad \text{and} \quad Y\_2 = h\_2(\mathcal{W}). \tag{7}$$

*The random variables Y*<sup>1</sup> *and Y*<sup>2</sup> *are counter-monotonic under* P*, if and only if,* (7) *holds with one of the functions h*1(·), *h*2(·) *being non-increasing, and the other non-decreasing.*

If two random variables are (counter) comonotonic under one probability measure, then they are also (counter) comonotonic under any other absolutely continuous probability measure, see, e.g., Proposition 2.1 of Cuestaalbertos et al. (1993). Thus, we omit the specification of the probability measure when discussing counter- and comonotonicity.

**Proposition 3** (Properties of Reverse Sensitivity Measure)**.** *The reverse sensitivity measure possesses the following properties:*


The function *s*(·) provides the flexibility to create sensitivity measures that quantify changes in moments, e.g., via *s*(*x*) = *x k* , *k* ∈ N, or in the tail of distributions, e.g., via *<sup>s</sup>*(*x*) = <sup>1</sup>{*x*>VaR*α*(*X<sup>i</sup>* )} , for *α* ∈ (0, 1).

Next, we generalise Definition 4 to a sensitivity measure that accounts for multiple input factors. While S Q∗ *<sup>i</sup>* measures the change of the distribution of *X<sup>i</sup>* from the baseline to the stressed model the sensitivity S Q∗ *i*,*j* introduced below, quantifies how the joint distribution of (*X<sup>i</sup>* , *Xj*) changes when moving from P to Q<sup>∗</sup> .

**Definition 6** (Bivariate Reverse Sensitivity Measure)**.** *For a function <sup>s</sup>*: <sup>R</sup><sup>2</sup> <sup>→</sup> <sup>R</sup>*, the reverse sensitivity measure to inputs* (*X<sup>i</sup>* , *Xj*) *with respect to a stressed probability measure* Q<sup>∗</sup> *is defined by*

$$\mathcal{S}\_{i,j}^{\mathbb{Q}^\*} = \begin{cases} \frac{\mathbb{E}^{\mathbb{Q}^\*} [s(\mathbf{X}\_i, \mathbf{X}\_j)] - \mathbb{E} [s(\mathbf{X}\_i, \mathbf{X}\_j)]}{\max\limits\_{\mathbb{Q} \in \mathcal{Q}} \mathbb{E}^{\mathbb{Q}} [s(\mathbf{X}\_i, \mathbf{X}\_j)] - \mathbb{E} [s(\mathbf{X}\_i, \mathbf{X}\_j)]} & \mathbb{E}^{\mathbb{Q}^\*} [s(\mathbf{X}\_i, \mathbf{X}\_j)] \ge \mathbb{E} [s(\mathbf{X}\_i, \mathbf{X}\_j)],\\ -\frac{\mathbb{E}^{\mathbb{Q}^\*} [s(\mathbf{X}\_i, \mathbf{X}\_j)] - \mathbb{E} [s(\mathbf{X}\_i, \mathbf{X}\_j)]}{\min\limits\_{\mathbb{Q} \in \mathcal{Q}} \mathbb{E}^{\mathbb{Q}} [s(\mathbf{X}\_i, \mathbf{X}\_j)] - \mathbb{E} [s(\mathbf{X}\_i, \mathbf{X}\_j)]} & \text{otherwise}, \end{cases}$$

*where* Q *is given in Definition 4.*

The bivariate sensitivity measure satisfies all the properties in Proposition 3 when *s*(*Xi*) is replaced by *s*(*X<sup>i</sup>* , *Xj*). The bivariate sensitivity S Q∗ *i*,*j* can also be generalised to *k* input factors by choosing a function *<sup>s</sup>*: <sup>R</sup>*<sup>k</sup>* <sup>→</sup> <sup>R</sup>.

**Remark 2.** *Probabilistic sensitivity measures are typically used for importance measurement and take values in* [0, 1]*; with 1 being the most important input factor and 0 being (desirably) independent from the output Borgonovo et al. (2021). This is in contrast to our framework where* S Q∗ *i lives in* [−1, 1] *and e.g., a negative dependence, such as negative quadrant dependence between s*(*Xi*) *and d*Q∗ *d*P *implies that* S Q∗ *<sup>i</sup>* < 0*, see Pesenti et al. (2019) [Proposition 4.3]. Thus, the proposed sensitivity measure is different in that it allows for negative sensitivities where the sign of* S Q∗ *i indicates the direction of the distributional change.*

#### **5. Application to a Spatial Model**

We consider a spatial model for modelling insurance portfolio losses where each individual loss occurs at different locations and the dependence between individual losses is a function of the distance between the locations of the losses. Mathematically, denote the locations of the insurance losses by *z*1, . . . , *z*10, where *z<sup>m</sup>* = (*z* 1 *<sup>m</sup>*, *z* 2 *<sup>m</sup>*) are the coordinates of location *zm*, *m* = 1, . . . , 10. The insurance loss at location *m*, denoted by *Lm*, follows a *Gamma*(5, 0.2 *m* ) distribution with location parameter 25. Thus, the minimum loss at each location is 25 and locations with larger mean also exhibit larger standard deviations. The losses *L*1, . . . , *L<sup>m</sup>* have, conditionally on Θ = *θ*, a Gaussian copula with correlation matrix given by *ρi*,*<sup>j</sup>* = Cor(*L<sup>i</sup>* , *Lj*) = *e* −*θ*||*zi*−*z<sup>j</sup>* ||, where || · || denotes the Euclidean distance. Thus, the further apart the locations *z<sup>i</sup>* and *z<sup>j</sup>* are the smaller the correlation between *L<sup>i</sup>* and *L<sup>j</sup>* . The parameter Θ takes values (0, 0.4, 5) with probabilities (0.05, 0.6, 0.35) that represent different regimes. Indeed, Θ = 0 corresponds to a correlation of 1 between all losses, independently of their location. Larger values of Θ correspond to smaller albeit still positive correlation. Thus, regime with Θ = 0 can be viewed as, e.g., circumstances suitable for natural disasters. We further define the total loss of the insurance company by *Y* = ∑ 10 *m*=1 *Lm*.

We perform two different stresses on the total loss *Y* detailed in Table 1. Specifically, we consider as a first stress a 0% change in HARA utility, a 0% change in ES0.8(*Y*), and a 1% increase in ES0.95(*Y*) from the baseline to the stressed model. The second stress is composed of a 1% increase in HARA utility, a 1% increase in ES0.8(*Y*), and a 3% increase in ES0.95(*Y*) compared to the baseline model. As the second stress increases all three metrics it may be viewed as a more severe distortion of the baseline model.


**Table 1.** Summary of the stresses applied to the portfolio loss *Y* represented in relative increases of the stressed model from the baseline model.

Next, we calculate reverse sensitivity measures for the losses *L*1, . . . , *L*<sup>10</sup> for both stresses Q∗ 1 and Q∗ 2 . Figure 5 displays the reverse sensitivity measures for functions *<sup>s</sup>*(*x*) = *<sup>x</sup>*, *<sup>s</sup>*(*x*) = <sup>1</sup>{*x*>*F*˘ *i* (0.8)} , and *<sup>s</sup>*(*x*) = <sup>1</sup>{*x*>*F*˘ *i* (0.95)} , from the left to the right panel, and where *F*˘ *i* , denotes the P-quantile function of *L<sup>i</sup>* , *i* = 1, . . . , 10.

**Figure 5.** Reverse sensitivity measures with *<sup>s</sup>*(*x*) = *<sup>x</sup>*, *<sup>s</sup>*(*x*) = <sup>1</sup>{*x*>*F*˘ *<sup>i</sup>*(0.8)} , and *<sup>s</sup>*(*x*) = <sup>1</sup>{*x*>*F*˘ *<sup>i</sup>*(0.95)} (left to right), for two different stresses on the output *Y*. First stress (salmon) is keeping the HARA utility and ES(*Y*)0.8 fixed and increasing the ES(*Y*)0.85 by 1%. Second stress (violet) is an increase of 1% in HARA utility, 1% ES(*Y*)0.8, and 3% in ES(*Y*)0.85.

We observe that for stress 2, the reverse sensitivities to all losses *L<sup>i</sup>* and all choices of function *s*(·) are positive. This contrasts the reverse sensitivities for stress 1. Indeed, for stress 1 the reverse sensitivities with both *<sup>s</sup>*(*x*) = *<sup>x</sup>* and *<sup>s</sup>*(*x*) = <sup>1</sup>{*x*>*F*˘ *i* (0.8)} are negative, with the former values being smaller indicating a smaller change in the distributions of the *Li* 's. By definition of the reverse sensitivity, the left panel corresponds to the (normalised) difference between the expectation under the stressed and baseline model. The middle and right panels correspond to the (normalised) change in the probability of exceeding *F*˘ *<sup>i</sup>*(0.8) and *F*˘ *<sup>i</sup>*(0.95), respectively. Thus, as seen in the plots, while the expectations and probabilities of exceeding the 80% P-quantile are smaller under the stressed model, the probabilities of exceeding the 95% P-quantile are increased substantially. The first stress increases the ES at level 0.95 while simultaneously fixes the utility and ES at level 0.8 to its values under the baseline model. This induces under the stressed probability measure a reduction of the mean and of the probability of exceeding the 80% P-quantile while the probability of exceeding the 95% P-quantile increases. Thus, the reverse sensitivity measures provide a spectrum of measures to analyse the distributional change of the losses *Li* from the baseline to the stressed model.

Next, for a comparison we calculate the delta sensitivity measure of introduced by Borgonovo (2007). For a probability measure Q the delta measure of *L<sup>i</sup>* is defined by

$$\mathcal{G}^{\mathbb{Q}}(L\_i) = \frac{1}{2} \int \int \left| f\_Y^{\mathbb{Q}}(y) - f\_{Y|i}^{\mathbb{Q}}(y \mid z) \right| f\_i^{\mathbb{Q}}(z) \, dy \, dz.$$

where *f* Q *Y* (·) and *f* Q *i* (·) are the densities of *Y* and *L<sup>i</sup>* under Q, respectively, and where *f* Q *Y*|*i* (·|·) is the conditional density of the total portfolio loss *Y* given *L<sup>i</sup>* under Q.

Table 2 reports the delta measures under the baseline model *ξ* <sup>P</sup> and the two stresses, i.e., *ξ* Q∗ <sup>1</sup> and *ξ* Q∗ <sup>2</sup> . We observe that the delta measures are similar for all losses *L<sup>i</sup>* and do not change significantly under the different probability measures. As the delta sensitivity measure quantifies the importance of input factors under a probability measure, having similar values for *ξ* P , *ξ* Q∗ <sup>1</sup> , and *ξ* Q∗ <sup>1</sup> , means that the importance ranking of the *L<sup>i</sup>* 's under different stresses does not change. We also report, in the first two columns of Table 2, the

reverse sensitivity measures with *<sup>s</sup>*(*x*) = <sup>1</sup>{*x*>*F*˘(0.95)} . The reverse sensitivity measures provide, in contrast to the delta measure, insight into the change in the distributions of the *Li* 's from P and Q∗ .


**Table 2.** Comparison of different sensitivity measures: First two columns correspond to the reverse sensitivity measures with *<sup>s</sup>*(*x*) = <sup>1</sup>{*x*>*F*˘(0.95)} and stressed models Q∗ 1 , and Q∗ 2 , respectively. The last three columns are the delta measure under P, Q∗ 1 , and Q∗ 2 , respectively.

Alternatively to considering the change in the marginal distributions *L<sup>i</sup>* from the baseline to the stressed model, we can study the change in the dependence between the losses when moving from the baseline to a stressed model. For this, we consider the bivariate reverse sensitivity measures for the pairs (*L*5, *L*10) and (*L*9, *L*10) for the second stress Q<sup>∗</sup> 2 , that is a 1% increase in HARA utility and ES0.8, and a 3% increase in ES0.95. Specifically, we look at the function *s*(*L<sup>i</sup>* , *<sup>L</sup>j*) = <sup>1</sup>{*Li*>*F*˘ *i* (0.95)} <sup>1</sup>{*Lj*>*F*˘ *j* (0.95)} , where *F*˘ *<sup>i</sup>*(·) and *<sup>F</sup>*˘ *<sup>j</sup>*(·) are the P-quantile functions of *L<sup>i</sup>* and *L<sup>j</sup>* respectively. This bivariate sensitivity measures quantifies the impact a stress has on the probability of joint exceedances with values S Q∗ 2 5,10 = 0.78 and S Q∗ 2 9,10 = 0.81 indicating that the probabilities of joint exceedances increase more for stress 2. This can also be seen in Figure 6 which shows the bivariate copulae contours of (*L*5, *L*10) (top panels) and (*L*9, *L*10) (bottom panels). The left contour plots correspond to the baseline model P whereas the right panels display the contours under the stress model Q∗ 2 (solid lines) together with the baseline contours (reported using partially transparent lines). The red dots are the simulated realisations of the losses (*L*5, *L*10) and (*L*9, *L*10), respectively (which are the same for the baseline and stressed model). We observe that for both pairs (*L*5, *L*10) and (*L*9, *L*10) the copula under the stressed model admits larger probabilities of joint large events, which is captured by the bivariate reverse sensitivity measure admitting positive values close to 1.

**Figure 6.** Contour plots of the bivariate copulae of (*L*5, *L*10) (**top** panels) and (*L*9, *L*10) (**bottom** panels) under different models. The left contour plots correspond to the baseline model and the right panels to the stress Q∗ 2 (solid lines) with the baseline contours reported using partially transparent lines. Red points are simulated realisations.

#### **6. Concluding Remarks**

We extend the reverse sensitivity analysis proposed by Pesenti et al. (2019) which proceeds as follows. Equipped with a baseline model which comprises of input and output random variables and a baseline probability measure, one derives a unique stressed model such that the output (or input) under the stressed model satisfies a prespecified stress and is closest to the baseline distribution. While Pesenti et al. (2019) consider the Kullback– Leibler divergence to measure the difference between the baseline and stressed models we utilise Wasserstein distance of order two. Compared to Pesenti et al. (2019) the Wasserstein distance allows for additional and different stresses on the output including the mean and variance, any distortion risk measure including the Value-at-Risk and Expected-Shortfall, and expected utility type constraints, thus making the reverse sensitivity analysis framework suitable for models used in financial and insurance risk management. We further discuss reverse sensitivity measures which quantify the change the inputs' distribution when moving from the baseline to a stressed model and illustrate our results on a spatial insurance portfolio application. The reverse sensitivity analysis framework (including the results from this work and from Pesenti et al. (2019) are implemented in the R package SWIM which is available on CRAN.

**Funding:** This research was funded by the Connaught Fund, the Canadian Statistical Sciences Institute (CANSSI), and the Natural Sciences and Engineering Research Council of Canada (NSERC) with funding reference numbers DGECR-2020-00333 and RGPIN-2020-04289.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** S.M.P. would like to thank Judy Mao for her help in implementing the numerical examples.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Proofs**

**Proof of Theorem 1.** We solve the optimisation on the set of quantile functions <sup>M</sup>˘ and define the Lagrangian with Lagrange multipliers *λ* = (*λ*1, . . . , *λ<sup>d</sup>* ) <sup>∈</sup> <sup>R</sup>*<sup>d</sup>*

$$\begin{split} L(\check{G}, \lambda) &= \int\_{0}^{1} \left( \check{G}(u) - \check{F}(u) \right)^{2} - 2 \sum\_{k=1}^{d} \lambda\_{k} \left( \check{G}(u) \gamma\_{k}(u) - r\_{k} \right) du \\ &= \int\_{0}^{1} \left( \check{G}(u) - \left( \check{F}(u) + \sum\_{k=1}^{d} \lambda\_{k} \gamma\_{k}(u) \right) \right)^{2} \\ &- 2 \sum\_{k=1}^{d} \lambda\_{k} \left( \check{F}(u) \gamma\_{k}(u) - r\_{k} \right) - \left( \sum\_{k=1}^{k} \lambda\_{k} \gamma\_{k}(u) \right)^{2} du. \end{split}$$

Thus, the optimisation problem (3) is equivalent to first solving, for fixed *λ*, the optimisation problem

$$\underset{\check{G}\in\check{\mathcal{M}}}{\text{arg min}} \, L(\check{G}, \lambda) \tag{A1}$$

and then finding *λ* such that the constraints are fulfilled. For fixed *λ*, the solution to (A1) is equal to the solution to

$$\underset{\check{G}\in\check{\mathcal{M}}}{\arg\min} \int\_{0}^{1} \left( \check{G}(u) - \left( \check{F}(u) + \sum\_{k=1}^{d} \lambda\_{k} \gamma\_{k}(u) \right) \right)^{2} du \,\,\mu$$

which is given by the isotonic projection of *F*˘(*u*) + ∑ *d k*=1 *λ<sup>k</sup> γ<sup>k</sup>* (*u*) onto the set <sup>M</sup>˘ and the Lagrange multipliers are such that the constraints are satisfied. Existence of the Lagrange multipliers follows since the set M is non-empty. Uniqueness follows by convexity of the Wasserstein distance, by convexity of the constraints on the set of quantile functions.

**Proof of Proposition 1.** For coherent distortion risk measures the corresponding weight function *γ* is non-decreasing. Moreover the optimal quantile function is given by Theorem 1 and is of the form *G*˘ *<sup>λ</sup>*(*u*) = *F*˘(*u*) + *λγ*(*u*) ↑ for some *λ* such that *G*˘ *<sup>λ</sup>* fulfils the constraint. The choice

$$
\lambda^\* = \frac{r - \rho\_\gamma(F)}{\int\_0^1 (\gamma(u))^2 du} \ge 0
$$

implies that *G*˘ *<sup>λ</sup>*<sup>∗</sup> (*u*) = *F*˘(*u*) + *λ* ∗*γ*(*u*) is a quantile function of the form (4) that fulfils the constraint. By uniqueness of Theorem 1, *G*˘ *<sup>λ</sup>*<sup>∗</sup> is indeed the unique solution.

**Proof of Theorem 2.** Since both constraints are convex in *G*˘ the Lagrangian with parameters *λ* = (*λ*1, . . . , *λ<sup>d</sup>* ) and *λ*˜ = (*λ*˜ <sup>1</sup>, . . . , *λ*˜ ˜*d* ) ≥ 0 becomes

$$\begin{split} L(\check{\mathbf{G}},\boldsymbol{\lambda},\check{\boldsymbol{\Lambda}}) &= \int\_{0}^{1} \left( \check{\mathbf{G}}(u) - \check{\mathbf{F}}(u) \right)^{2} + 2 \sum\_{k=1}^{d} \lambda\_{k} \left( h\_{k}(u) \check{\mathbf{G}}(u) - c\_{k} \right) \, du \\ &+ \sum\_{k=1}^{\tilde{d}} \tilde{\lambda}\_{k} \Big( \tilde{h}\_{k}(u) \left( \check{\mathbf{G}}(u) \right)^{2} - \check{\mathbf{c}}\_{k} \Big) \, du \\ &= \int\_{0}^{1} \tilde{\Lambda}(u) \left( \check{\mathbf{G}}(u) - \frac{1}{\check{\mathbf{A}}(u)} \left( \check{\mathbf{F}}(u) - \sum\_{k=1}^{d} \lambda\_{k} h\_{k}(u) \right) \right)^{2} \\ &- \frac{1}{\check{\mathbf{A}}(u)} \left( \check{\mathbf{F}}(u) - \sum\_{k=1}^{d} \lambda\_{k} h\_{k}(u) \right)^{2} + \left( \check{\mathbf{F}}(u) \right)^{2} - 2 \sum\_{k=1}^{d} \lambda\_{k} \boldsymbol{c}\_{k} - \sum\_{k=1}^{d} \check{\lambda}\_{k} \check{\mathbf{c}}\_{k} \, du \end{split}$$

where Λ˜ (*u*) = 1 + ∑ ˜*d k*=1 *λ*˜ *k* ˜*hk* (*u*). Since *<sup>λ</sup>*˜ <sup>≥</sup> <sup>0</sup> by the KKT-condition, we obtain that <sup>Λ</sup>˜ (*u*) <sup>≥</sup> <sup>0</sup> for all *<sup>u</sup>* <sup>∈</sup> [0, 1]. Therefore, for fixed *<sup>λ</sup>*, *<sup>λ</sup>*˜ , using an argument similar to the proof of Theorem 1, the solution (as a function of *λ*, *λ*˜ ) is given by the weighted isotonic projection of <sup>1</sup> Λ˜ (*u*) *<sup>F</sup>*˘(*u*) <sup>−</sup> <sup>∑</sup> *d k*=1 *λkh<sup>k</sup>* (*u*) , with weight function <sup>Λ</sup>˜ (·).

**Proof of Proposition 2.** The mean and variance constraint can be rewritten as

$$\begin{aligned} m' &= \int \mathbf{x} \, d\mathbf{G}(\mathbf{x}) = \int\_0^1 \breve{\mathbf{G}}(\boldsymbol{u}) \, d\boldsymbol{u} \qquad \text{and} \\ \left(\boldsymbol{\sigma'}\right)^2 &= \int \left(\mathbf{x} - \boldsymbol{m'}\right)^2 \, d\mathbf{G}(\mathbf{x}) = \int\_0^1 \left(\breve{\mathbf{G}}(\boldsymbol{u}) - \boldsymbol{m'}\right)^2 \, d\boldsymbol{u} \, d\boldsymbol{u} \end{aligned}$$

Thus, the Lagrangian with Lagrange multipliers *λ* = (*λ*1, . . . , *λk*+<sup>2</sup> ) is, if *λ*<sup>2</sup> 6= −1,

$$\begin{split} L(\check{\mathbf{G}},\lambda) &= \int\_{0}^{1} \left( \check{\mathbf{G}}(u) - \check{\mathbf{F}}(u) \right)^{2} du - 2\lambda\_{1} \left( \int\_{0}^{1} \check{\mathbf{G}}(u) \, du - m' \right) \\ &\quad + \lambda\_{2} \Big( \int\_{0}^{1} \left( \check{\mathbf{G}}(u) - m' \right)^{2} du - \left( \sigma' \right)^{2} \Big) \\ &\quad - 2 \sum\_{k=1}^{d} \lambda\_{k+2} \Big( \int\_{0}^{1} \check{\mathbf{G}}(u) \gamma\_{\mathbf{K}}(u) \, du - r\_{k} \Big) \\ &= \left( 1 + \lambda\_{2} \right) \int\_{0}^{1} \left( \check{\mathbf{G}}(u) - \frac{1}{1 + \lambda\_{2}} \left( \check{\mathbf{F}}(u) + \lambda\_{1} + \lambda\_{2} m' + \sum\_{k=1}^{d} \lambda\_{k+2} \gamma\_{\mathbf{K}}(u) \right) \right)^{2} \\ &\quad - \frac{1}{1 + \lambda\_{2}} \left( \check{\mathbf{F}}(u) + \lambda\_{1} + \lambda\_{2} m' + \sum\_{k=1}^{d} \lambda\_{k+2} \gamma\_{\mathbf{K}}(u) \right)^{2} + \left( \check{\mathbf{F}}(u) \right)^{2} \, du \\ &\quad + 2\lambda\_{1} m' + \lambda\_{2} \Big( \left( m' \right)^{2} - \left( \sigma' \right)^{2} \Big) + 2 \sum\_{k=1}^{d} \lambda\_{k+2} \tau\_{k} .\end{split}$$

For fixed Lagrange multipliers *λ* with *λ*<sup>2</sup> 6= −1, the optimal quantile function is characterised by the isotonic projection and given by (using an analogous argument to the proof of Theorem 1)

$$\begin{split} \check{G}^\*(u) &= \left( \frac{1}{1+\lambda\_2} \Bigl( \check{F}(u) + \lambda\_1 + \lambda\_2 m' + \sum\_{k=1}^d \lambda\_{k+2} \gamma\_k(u) \Bigr) \right)^\uparrow \\ &= \frac{1}{|1+\lambda\_2|} \Bigl( \text{sgn}(1+\lambda\_2) \Bigl( \check{F}(u) + \lambda\_1 + \lambda\_2 m' + \sum\_{k=1}^d \lambda\_{k+2} \gamma\_k(u) \Bigr) \Bigr)^\uparrow \\ &= \frac{1}{|1+\lambda\_2|} \Bigl( \check{H}(u) \Bigr) . \end{split} \tag{A2}$$

where we define *<sup>H</sup>*˘ (*u*) = sgn(1 + *λ*2) *F*˘(*u*) + *λ*<sup>1</sup> + *λ*2*m*<sup>0</sup> + ∑ *d k*=1 *λk*+<sup>2</sup> *γ<sup>k</sup>* (*u*) ↑ <sup>∈</sup> <sup>M</sup>˘ , and sgn(·) denotes the sign function. Next we show that *λ*<sup>2</sup> cannot be in a neighbourhood of −1. It holds that for *λ*<sup>2</sup> 6= −1,

$$\int\_0^1 \left(\check{\mathcal{G}}^\*(u)\right)^2 du = \frac{1}{(1+\lambda\_2)^2} \int\_0^1 \left(\check{H}(u)\right)^2 du\,. \tag{A3}$$

Since the rhs of (A3) is increasing for | *λ*<sup>2</sup> + 1 | & 0, there exists a *ε*<sup>0</sup> > 0 such that for all *ε* < *ε*<sup>0</sup> and *λ*<sup>2</sup> ∈ (−1 − *ε*, −1 + *ε*), it holds that

$$\frac{1}{(1+\lambda\_2)^2} \int\_0^1 \left(\tilde{H}(u)\right)^2 du \,> \left(\sigma'\right)^2 + \left(m'\right)^2 \lambda\_2$$

which is a contradiction to the optimality of *G*˘ <sup>∗</sup> . Thus, *λ*<sup>2</sup> is indeed bounded away from −1 and the unique solution is given in (A2).

**Proof of Theorem 3.** We split this proof into the two cases (*i*), that is constraint (a) and (*ii*), i.e., constraint (b).

Case (*i*): For constraint a), i.e., VaR*α*(*G*) = *q*, we first assume that *q* ≤ VaR*α*(*F*) which implies *<sup>F</sup>*˘(*αF*) = *<sup>q</sup>* <sup>≤</sup> *<sup>F</sup>*˘(*α*) and thus *<sup>α</sup><sup>F</sup>* <sup>≤</sup> *<sup>α</sup>*. Therefore, *<sup>G</sup>*˘ <sup>∗</sup> (*u*) = *F*˘(*u*) + *q* − *F*˘(*u*) <sup>1</sup>{*u*∈(*αF*,*α*]} is a quantile function which satisfies the constraint. Next, we show that *G* ∗ has a smaller Wasserstein distance to *F* than any other distribution function satisfying the constraint. For this, let *<sup>H</sup>*˘ be a quantile function satisfying the constraint and *<sup>H</sup>*˘ (*u*) <sup>6</sup><sup>=</sup> *<sup>G</sup>*˘(*u*) on a measurable set of non-zero measure. Then

$$\begin{aligned} \mathcal{W}\_2(H, F) &= \int\_0^{a\_F} \left( \check{\mathcal{H}}(u) - \check{\mathcal{F}}(u) \right)^2 du + \int\_{a\_F}^a \left( \check{\mathcal{H}}(u) - \check{\mathcal{F}}(u) \right)^2 du + \int\_a^1 \left( \check{\mathcal{H}}(u) - \check{\mathcal{F}}(u) \right)^2 du \\ &\ge \int\_{a\_F}^a \left( \check{\mathcal{H}}(u) - \check{\mathcal{F}}(u) \right)^2 du \,. \end{aligned}$$

By non-decreasingness of *<sup>H</sup>*˘ and *<sup>F</sup>*˘ and by the constraint it holds for all *<sup>u</sup>* <sup>∈</sup> [*αF*, *<sup>α</sup>*] that *<sup>H</sup>*˘ (*u*) <sup>≤</sup> *<sup>H</sup>*˘ (*α*) = *<sup>q</sup>* <sup>=</sup> *<sup>F</sup>*˘(*αF*) <sup>≤</sup> *<sup>F</sup>*˘(*u*). Thus, on the interval [*αF*, *<sup>α</sup>*], we obtain *<sup>H</sup>*˘ (*u*) <sup>−</sup> *F*˘(*u*) <sup>2</sup> <sup>≥</sup> *<sup>q</sup>* <sup>−</sup> *<sup>F</sup>*˘(*u*) 2 and therefore

$$\mathcal{W}\_2(H, \mathbb{F}) \ge \int\_{a\_F}^{a} \left(\check{H}(u) - \check{\mathbb{F}}(u)\right)^2 du \ge \int\_{a\_F}^{a} \left(q - \check{\mathbb{F}}(u)\right)^2 du = \mathcal{W}\_2(G^\*, \mathbb{F}) \ge 0$$

where at least one inequality is strict since *<sup>H</sup>*˘ (*u*) <sup>6</sup><sup>=</sup> *<sup>G</sup>*˘(*u*) on a measurable set of non-zero measure. Uniqueness follows by the strict convexity of the Wasserstein distance and since the constraint is convex on the set of quantile functions.

Second, we assume that *q* > VaR*α*(*F*) and show that there does not exist a solution. Assume by contradiction that *G*˘ is an optimal quantile function satisfying the constraint. By definition of *<sup>α</sup>F*, we have that *<sup>q</sup>* <sup>=</sup> *<sup>F</sup>*˘(*αF*) <sup>&</sup>gt; *<sup>F</sup>*˘(*α*) and thus *<sup>α</sup><sup>F</sup>* <sup>≥</sup> *<sup>α</sup>*. We apply a similar argument to the first part of the proof using non-decreasingness of *G*˘, *G*˘(*α*) = *q*, and optimality of *G*˘, to obtain that *G*˘ is constant equal to *q* on [*α*, *αF*] and equal to *F*˘(*u*) for *u* > *αF*. Specifically, it holds that

$$\check{G}(\mu) = \check{F}(\mu) + (q - \check{F}(\mu))\mathbb{1}\_{\{\mu \in (a, \alpha\_F]\}\prime} \quad \text{for all} \quad \mu > \alpha\text{ .}$$

Moreover, since the optimal quantile function minimises the Wasserstein distance to *F*, it holds that, for all *ε* > 0, *G*˘ satisfies

$$
\vec{G}(\mathfrak{u}) = \vec{F}(\mathfrak{u})\,, \quad \text{for all} \quad \mathfrak{u} \le \mathfrak{a} - \varepsilon \; .
$$

Thus, we can define for all *ε* ∈ (0, *α*) the family of quantile functions

$$\check{H}\_{\varepsilon}(\mu) = \check{F}(\mu) + (q - \check{F}(\mu))\mathbb{1}\_{\{\mu \in (\alpha - \varepsilon, \kappa\_F]\}}\star$$

which satisfies *W*2(*Hε*<sup>1</sup> , *F*) < *W*2(*Hε*<sup>2</sup> , *<sup>F</sup>*) for all <sup>0</sup> <sup>≤</sup> *<sup>ε</sup>*<sup>1</sup> <sup>&</sup>lt; *<sup>ε</sup>*2, and *<sup>H</sup>*˘ *<sup>ε</sup>*(*α*) = *q* for all *ε* > 0. However, lim*ε*&<sup>0</sup> *<sup>H</sup>*˘ *<sup>ε</sup>*(*α*) = *<sup>F</sup>*˘(*α*) <sup>&</sup>lt; *<sup>q</sup>* and thus the quantile function lim*ε*&<sup>0</sup> *<sup>H</sup>*˘ *<sup>ε</sup>*(*u*) does not fulfil the constraint. Hence, we obtain a contradiction to the optimality of *G*˘.

Case (*ii*): First, we assume that *<sup>q</sup>* <sup>≥</sup> VaR<sup>+</sup> *α* (*F*) which implies that *<sup>F</sup>*˘(*αF*) = *<sup>q</sup>* <sup>≥</sup> *<sup>F</sup>*˘+(*α*) <sup>≥</sup> *<sup>F</sup>*˘(*α*) and thus *<sup>α</sup><sup>F</sup>* <sup>≥</sup> *<sup>α</sup>*. Therefore *<sup>G</sup>*˘ <sup>∗</sup> (*u*) = *F*˘(*u*) + *<sup>q</sup>* <sup>−</sup> *<sup>F</sup>*˘(*u*) <sup>1</sup>{*u*∈(*α*,*αF*]} is a quantile function. Moreover, *G*˘ <sup>∗</sup> satisfies the constraint since by right-continuity of *G*˘ <sup>∗</sup> , we have that

$$\check{G}^{\*+}(\mathfrak{a}) = \lim\_{\varepsilon \searrow 0} \check{G}^{\*+}(\mathfrak{a} + \varepsilon) = q \cdot \mathfrak{a}$$

The proof that *G*˘ <sup>∗</sup> has the smallest Wasserstein distance to *F* compared to any other distribution function satisfying the constraint is analogous to the one in case (*i*).

For the case when *q* > VaR<sup>+</sup> *α* (*F*), the argument of non-existence of the solution follows using similar arguments as those in case (*i*).

**Proof of Theorem 4.** By concavity of the utility function, the constraint is convex and can be written as − R 1 0 *u G*˘(*v*) *dv* + *c* ≤ 0. Thus, we can define the Lagrangian with *λ*<sup>1</sup> ≥ 0 and (*λ*2, . . . , *λd*+<sup>1</sup> ) <sup>∈</sup> <sup>R</sup>*<sup>d</sup>* by

$$\begin{aligned} L(\check{G}, \lambda) &= \frac{1}{2} \int\_0^1 \left( \check{G}(v) - \check{F}(v) \right)^2 - \lambda\_1 \left( u \left( \check{G}(v) \right) - c \right) - \sum\_{k=1}^d \lambda\_{k+1} \left( \check{G}(v) \gamma\_k(v) - r\_k \right) dv \\ &= \int\_0^1 T(\check{G}(v)) - \check{G}(v) \left( \check{F}(v) + \sum\_{k=1}^d \lambda\_{k+1} \gamma\_k(v) \right) \\ &+ \frac{1}{2} \left( \check{F}(v) \right)^2 + \lambda\_1 c + \sum\_{k=1}^d \lambda\_{k+1} r\_k \, dv \end{aligned}$$

where *T*(*x*) = <sup>1</sup> 2 *x* <sup>2</sup> <sup>−</sup> *<sup>λ</sup>*1*u*(*x*). Therefore, for fixed *<sup>λ</sup>*1, . . . , *<sup>λ</sup>d*+<sup>1</sup> , we apply Theorem 3.1 by Barlow and Brunk (1972) and obtain the unique optimal quantile function (as a function of *λ*1, . . . , *λd*+<sup>1</sup> ), that is *G*˘ <sup>∗</sup> (*v*) = *ν*˘*λ*<sup>1</sup> *F*˘(*u*) + ∑ *d k*=1 *λk*+1*γ<sup>k</sup>* (*v*) ↑ , where *ν*˘*λ*<sup>1</sup> is the left-inverse of *νλ*<sup>1</sup> (*x*) = *x* − *λ*<sup>1</sup> *u* 0 (*x*).

Next, we show that if *d* = 0, the utility constraint is binding, that is *λ*<sup>1</sup> > 0. For this, assume by contradiction that the *λ*<sup>1</sup> = 0, then the optimal quantile function becomes *G*˘ ∗ (*u*) = *ν*˘<sup>0</sup> *F*˘(*u*) . Since *ν*0(*x*) = *x*, we obtain that *G*˘ <sup>∗</sup> (*u*) = *F*˘(*u*). *F*˘, however, does not fulfil the constraint, which is a contradiction to the optimality of *G*˘ <sup>∗</sup> .

**Proof of Proposition 3.** We prove the properties one-by-one:

(*i*) We first define for a random variable *Z* with P-distribution *F<sup>Z</sup>* the random variable *U<sup>Z</sup>* := *FZ*(*Z*). Then, *U<sup>Z</sup>* and *Z* are comonotonic and *U<sup>Z</sup>* has a uniform distribution under P. Next, recall that for any random variables *Y*1,*Y*<sup>2</sup> it holds that Rüschendorf (1983)

$$\mathbb{E}\left[\mathbf{Y}\_1\,\mathbf{F}\_{\mathbf{Y}\_2}^{-1}\left(1-\mathcal{U}\_{\mathbf{Y}\_1}\right)\right] \le \mathbb{E}[\mathbf{Y}\_1\,\mathbf{Y}\_2] \le \mathbb{E}\left[\mathbf{Y}\_1\,\mathbf{F}\_{\mathbf{Y}\_2}^{-1}\left(\mathcal{U}\_{\mathbf{Y}\_1}\right)\right].\tag{A4}$$

where *F* −1 *Y*2 *UY*<sup>1</sup> is the random variable that is comonotonic to *Y*<sup>1</sup> and has the same P-distribution as *Y*2. Similarly, *F* −1 *Y*2 1 − *UY*<sup>1</sup> is the random variable that is countermonotonic to *Y*<sup>1</sup> and has the same P-distribution as *Y*2. The left (right) inequality in (A4) become equality if and only if the random variables *Y*<sup>1</sup> and *Y*<sup>2</sup> are countercomonotonic (comonotonic).

Thus, we can rewrite the maximum in the normalising constant of the reverse sensitivity measure as follows

$$\max\_{\mathbb{Q}\in\mathcal{Q}}\mathbb{E}^{\mathbb{Q}}[s(X)] = \max\_{\mathbf{Z}\stackrel{\mathbb{P}}{=}\frac{d\mathbb{Q}^\*}{d\mathbb{P}}}\mathbb{E}[s(X)\,\mathbf{Z}] = \mathbb{E}\left[s(X)\,\mathcal{F}\_{\frac{d\mathbb{Q}^\*}{d\mathbb{P}}}^{-1}\left(\mathcal{U}\_{s(X)}\right)\right].$$

and the minimum in the normalising constant is

$$\min\_{\mathcal{Q}\in\mathcal{Q}}\mathbb{E}^{\mathcal{Q}}[s(X)] = \min\_{Z\stackrel{\mathbb{P}}{=}\frac{d\mathcal{Q}^\*}{d\mathbb{P}}}\mathbb{E}[s(X)\,Z] = \mathbb{E}\left[s(X)\,F\_{\frac{d\mathcal{Q}^\*}{d\mathbb{P}}}^{-1}\left(1-\mathcal{U}\_{s(X)}\right)\right].$$

The reverse sensitivity for the case E Q∗ [*s*(*Xi*)] ≥ E[*s*(*Xi*)] then becomes

$$\mathcal{S}\_i^{\mathbb{Q}^\*} = \frac{\mathbb{E}[s(X\_i)\frac{d\mathbb{Q}^\*}{d\mathbb{P}}] - \mathbb{E}[s(X\_i)]}{\mathbb{E}\left[s(X)\,F\_{\frac{d\mathbb{Q}^\*}{d\mathbb{P}}}\left(\mathcal{U}\_{s(X)}\right)\right] - \mathbb{E}[s(X\_i)]} \,\,^\*$$

which satisfies 0 ≤ *S* Q∗ *<sup>i</sup>* ≤ 1 using again (A4). For the case E Q∗ [*s*(*Xi*)] ≤ E[*s*(*Xi*)], it holds that

$$S\_i^{\mathbb{Q}^\*} = -\frac{\mathbb{E}[s(X\_i)\frac{d\mathbb{Q}^\*}{d\mathbb{P}}] - \mathbb{E}[s(X\_i)]}{\mathbb{E}\left[s(X)F\_{\frac{d\mathbb{Q}^\*}{d\mathbb{P}}}\left(1 - \mathcal{U}\_{s(X)}\right)\right] - \mathbb{E}[s(X\_i)]} \geqslant$$

which satisfies −1 ≤ *S* Q∗ *<sup>i</sup>* ≤ 0.

(*ii*) Assume that *s*(*Xi*) and *<sup>d</sup>*Q<sup>∗</sup> *d*P are independent under P, then

 $\mathbb{E}[s(X\_i)\frac{d\mathbb{Q}^\*}{d\mathbb{P}}] = \mathbb{E}[s(X\_i)]$  $\mathbb{E}\left[\frac{d\mathbb{Q}^\*}{d\mathbb{P}}\right] = \mathbb{E}[s(X\_i)]$ 

and the reverse sensitivity measure is indeed zero.


The proof that the joint reverse sensitivity *S* Q∗ *i*,*j* also fulfils the above properties follows using analogous arguments and replacing *s*(*Xi*) with *s*(*X<sup>i</sup>* , *Xj*).

#### **References**

Acerbi, Carlo, and Dirk Tasche. 2002. On the coherence of Expected Shortfall. *Journal of Banking & Finance* 26: 1487–503.


Barlow, Richard E., David J. Bartholomew, John. M. Bremner, and Hugh D. Brunk. 1972. *Statistical Inference under Order Restrictions: The Theory and Application of Isotonic Regression*. Hoboken: Wiley.

Bernard, Carole, Silvana M. Pesenti, and Steven Vanduffel. 2020. Robust distortion risk measures. *arXiv* arXiv:2205.08850.


Rüschendorf, Ludger. 1983. Solution of a statistical optimization problem by rearrangement methods. *Metrika* 30: 55–61. [CrossRef] Saltelli, Andrea, Marco Ratto, Terry Andres, Francesca Campolongo, Jessica Cariboni, Debora Gatelli, Michaela Saisana, and Stefano Tarantola. 2008. *Global Sensitivity Analysis: The Primer*. Hoboken: John Wiley & Sons.


## *Article* **Special-Rate Life Annuities: Analysis of Portfolio Risk Profiles**

**Ermanno Pitacco \* and Daniela Y. Tabakova**

Insurance and Risk Management, MIB Trieste School of Management, Largo Caduti di Nasiriya 1, 34142 Trieste, Italy; tabakova@mib.edu

**\*** Correspondence: ermanno.pitacco@deams.units.it

**Abstract:** Special-rate life annuities are life annuity products whose single premium is based on a mortality assumption driven (at least to some extent) by the health status of the applicant. The health status is ascertained via an appropriate underwriting step (which explains the alternative expression "underwritten life annuities"). Better annuity rates are then applied in presence of poor health conditions. The worse the health conditions, the smaller the modal age at death (as well as the expected lifetime), but the higher the variance of the lifetime distribution. The latter aspect is due to significant data scarcity as well as to the mix of possible pathologies leading to each specific rating class. A higher degree of (partially unobservable) heterogeneity inside each sub-portfolio of special-rate annuities follows, and this results in a higher variability of the total portfolio payout. The present research aims at analyzing the impact of extending the life annuity portfolio by selling special-rate life annuities. Numerical evaluations have been performed by adopting a deterministic approach as well as a stochastic one, according to diverse assumptions concerning both lifetime distributions and portfolio structure and size. Our achievements witness the possibility of extending the annuity business without taking huge amounts of risk. Hence, the risk management objective "enhancing the company market share" can be pursued without significant worsening of the annuity portfolio risk profile.

**Keywords:** life annuities; standard annuities; underwritten annuities; enhanced annuities; impaired annuities; preferred risks; substandard lives

#### **1. Introduction and Motivation**

Considerable attention is currently being devoted in insurance work (and, in particular, in the actuarial work) to the management of life annuity portfolios and to the annuity product design, because of the growing importance of annuity benefits paid by private pension schemes and individual policies.

In particular, the progressive shift in many countries from defined benefit to defined contribution pension schemes has increased the interest in life annuity products with a guaranteed periodic benefit. Nevertheless, various "weak" features of the (standard) life annuities should be noted, looking at the product from both the annuity provider's and the customer's perspective.

However, many features can be improved by moving from traditional products to more complex products, for example, by adding riders (that is, supplementary benefits), or by adopting restrictions on the age intervals covered, or by allowing for individual risk factors; hence, "tailoring" the annuity rates (at least to some extent) to specific features of the customer.

Special-rate life annuities are life annuity products whose single premium is based on a mortality assumption driven by the health status of the applicant. The health status is ascertained via an appropriate underwriting step (which explains the alternative expression "underwritten life annuities"). Better annuity rates are then applied in presence of poor health conditions. The worse the health conditions, the smaller the modal age at death (as well as the expected lifetime), but the higher the variance of the lifetime distribution. The

**Citation:** Pitacco, Ermanno, and Daniela Y. Tabakova. 2022. Special-Rate Life Annuities: Analysis of Portfolio Risk Profiles. *Risks* 10: 65. https://doi.org/10.3390/ risks10030065

Academic Editor: Nadine Gatzert

Received: 5 February 2022 Accepted: 10 March 2022 Published: 13 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

latter aspect is due to significant data scarcity as well as to the mix of possible pathologies leading to each specific rating class. A higher degree of (partially unobservable) heterogeneity inside each sub-portfolio of special-rate annuities follows, and this results in a higher variability of the total portfolio payout. By selling special-rate annuities, on the one hand, a higher premium income can be expected, and on the other, a higher variability of the portfolio payout must be faced. What about the "balance"? Our achievements witness the possibility of extending the annuity business without taking huge amounts of risk. Hence, the risk management objective "enhancing the company market share" can be pursued without significant worsening of the annuity portfolio risk profile.

Diverse input data might lead to worse risk profiles. An appropriate sensitivity testing can then help in checking risk profile changes. Assuming the extension of the life annuity business as the insurer's target, the present paper aims at providing a simple technical tool for assessing how and to what extent selling special-rate annuities impacts the portfolio risk profile. The structure of the proposed tool arises from a trade-off between strictly pragmatic approaches (frequently adopted in current actuarial practice) and rigorous mathematical settings (which may result in implementation difficulties, notably because of data scarcity).

The remainder of the paper is organized as follows. A compact literature review is provided in Section 2, while Section 3 describes the main products in the area of special-rate annuities. Biometric assumptions underlying the assessment of portfolio risk profiles are defined and commented in Section 4. In Section 5, portfolio structures are specified in terms of sub-portfolio sizes; then, results of interest, which can express the portfolio risk profile, are defined. Our main achievements are discussed in Sections 6 and 7 where numerical results obtained by adopting a deterministic and a stochastic approach are respectively presented. Finally, Section 8 concludes the paper.

#### **2. Literature Review**

The main features of life annuity products are discussed in many life insurance and actuarial textbooks. A presentation of basic actuarial models for premium and reserve calculations for life annuities as well as a discussion of possible innovations in life annuity products are provided by Pitacco (2021).

Heterogeneity in mortality and risk classification constitute the natural frameworks in which the basic features of special-rate life annuities can be analyzed. Risk classification in life insurance and life annuities is addressed in many books and papers; a compact review, together with an extensive reference list, is provided by Haberman and Olivieri (2014). The impact of risk classification on the structure of life annuity portfolios is dealt with by Gatzert et al. (2012), Hoermann and Russ (2008) and Olivieri and Pitacco (2016).

The impact of heterogeneity on portfolio results and the consequent capital requirements are analyzed by Denuit and Frostig (2006). More specifically, Denuit and Frostig (2007) focusses on heterogeneity among lifetimes in the context of stochastic mortality according to the Lee-Carter model.

Heterogeneity in mortality is due to both observable and unobservable risk factors. The reader can refer to Pitacco (2019) for a literature review from an actuarial perspective, as well as for a discussion of models, which can be used to represent specific mortality rates accounting for observable risk factors.

Special-rate life annuities are described in various papers and technical reports: see, in particular Ainslie (2000), Drinkwater et al. (2006), Ridsdale (2012) and Rinke (2002). The article by Edwards (2008) is specifically devoted to life annuity rating based on the postcode (that is, a proxy for social class and location of housing).

An interesting analysis of market issues related to special-rate annuities is presented by Gatzert and Klotzki (2016), where barriers on the supply side and the demand side are in particular addressed. Practical aspects of pricing special-rate life annuities are dealt with by Gracie and Makin (2006) and James (2016).

Special-rate annuities in the context of new product development (NPD) processes are addressed in Chapter 9 of Pitacco (2020). In particular, the NPD according to the structure of the risk management process and the logic of the Stage-Gate® process is described and commented on<sup>1</sup> .

Underwriting for special-rate life annuities can be implemented in a number of ways, and several classifications can be conceived. An interesting classification has been proposed by Rinke (2002), and summarized in Chapter 9 of Pitacco (2021).

Statistics regarding extra-mortality by various causes are beyond the scope of this paper. Here, we only cite the contribution by Weinert (2006), which has suggested some baseline choices for the lifetime distributions. For a list of references the reader can refer to Pitacco and Tabakova (2020).

#### **3. The Products**

The terminology adopted in the technical literature to denote the various types of special-rate annuities is not univocally defined. For example, the term "enhanced annuity" is frequently used in the wider sense of all life annuities where, given the single premium, the annual benefit is of a higher level than the standard, due to some customer's characteristics (health, lifestyle, etc.). In what follows, we refer to the terminology originally adopted to define special-rate annuities, which are sold in the UK market (see Ridsdale 2012). The same terminology has also been adopted in the book (Pitacco 2021), from which the following descriptions have been taken.

The following special-rate annuities are sold in several markets.

	- (a) *Smoker life annuities*: if the applicant has smoked at least a given number of cigarettes for a certain number of years, then they are eligible for a smoker annuity.
	- (b) Mortality differences between married and unmarried individuals underpin the use of special rates in pricing the *unmarried lives annuities*. The observed higher mortality rates of unmarried individuals justify a higher annuity rate.

Thus, moving from type 1 to type 4 results in progressively higher mortality assumptions, shorter life expectancy, and hence, for a given single premium amount, in higher annuity benefits. Of course, an insurer can decide to offer a more limited set of products.

The applicant's health status and, notably, the presence of past or current diseases is explicitly considered in the special-rate annuities of types 2, 3 and 4. Various factors concerning the health status can be accounted for, and medical ascertainment is of course required. In particular, the underwriting process for impaired-life annuities and care annuities must result in classifying the applicant as a *substandard risk*, because of ascertainment of significant extra-mortality. For this reason, annuities of types 3 and 4 are sometimes named *substandard life annuities*.

The above list of special-rate annuity types can be completed by the *postcode life annuities*, which constitute an interesting example of "environment-based" rating. The postcode can provide a proxy for social class and location of housing; that is, risk factors that may have a significant impact on the lifestyle and hence on the life expectancy. Then, its use as a rating factor for pricing life annuities can be justified.

#### **4. The Mortality Model**

To assess present values (that is, random present values and expected present values) of benefits paid by special-rate annuities, a model quantifying the mortality of annuitants (standard annuitants and special-rate annuitants) must be chosen. The following aspects must be considered in particular:


#### *4.1. General Aspects*

A higher degree of heterogeneity in mortality affects a life annuity portfolio also including special-rate annuities, with respect to a standard annuity portfolio. As is well known, heterogeneity may be due to both observable and unobservable risk factors, which imply diverse modeling choices. In the case of a life annuity portfolio, the practical problem is: to what extent the underwriting process can detect, for each applicant, the outcomes of the risk factors that entitle the individual to purchase a special-rate annuity? Even a rigorous underwriting process can leave some degree of "residual" heterogenity inside each special-rate class, because:


Given that some degree of heterogeneity inside each risk class is unavoidable, the problem is how to model its impact on the individual age-pattern of mortality. We focus on the two following choices.

The modeling of heterogeneity in mortality due to unobservabke risk factors has found an elegant and rigorous solution in the concept of (constant) frailty, initially described by Beard (1959), but formally defined by Vaupel et al. (1979). A well known implementation of frailty modeling leads to the so-called Gompertz–Gamma model, which results in one of the Perks laws (see Perks 1932). A number of generalizations of the concept of frailty have been proposed, in particular looking at possible dynamic features of the individual frailty (a survey is presented, for example, by Pitacco 2019).

Frailty modeling to express heterogeneity in mortality has been adopted by various Authors: see, for example, Denuit and Frostig (2007), and, in the context of special-rate annuities, by Olivieri and Pitacco (2016).

Given the practical difficulties in calibrating the frailty model, the paper by Olivieri and Pitacco (2016) specifically aims at assessing, via sensitivity analysis, the impact of diverse assumptions for the frailty parameter values (notably, the Gamma parameters) on portfolio results of interest.

A simpler choice can conversely consist in implicitly expressing the presence of heterogeneity in mortality by directly assuming a higher variance in the individual lifetime distribution, as suggested by statistical data analyzed by various Authors (see, for example, Weinert 2006), that is, the stronger the assumed degree of heterogeneity, the higher the variance. If a mortality law is chosen to express the individual age-pattern of mortality, a sensitivity analysis can be performed also in this setting.

As far as the future mortality trends are concerned, the following aspects can be singled out:


In what follows, we express heterogeneity inside each rating class via the variance of the individual lifetime distribution (see Section 4.2). We disregard systematic longevity risk, as heterogeneity mainly affects the idiosyncratic risk in each rating class; that is, the risk of random fluctuations around the expected value.

Assumptions regarding the relations among individual lifetimes must supplement the mortality model. Given the above hypotheses, in what follows we assume that individual lifetimes are independent random variables. While this is, of course, a simplifying assumption, correlation among lifetimes could be assumed, in particular to express uncertainty about future mortality trends in the context of a stochastic mortality model. Conditional independence would replace, in that case, the independence assumption.

#### *4.2. Age Pattern of Mortality*

Three (hypothetical) curves showing expected number *d<sup>x</sup>* of deaths between exact age *x* and *x* + 1 out of a notional cohort of 100,000 individuals are shown in Figure 1. We recall that the expected numbers of deaths are proportional to the values of the probability density function in the time-continuous context defined below.

**Figure 1.** Curves of deaths for different life annuities.

The worse the health conditions, the smaller the modal age at death (as well as the life expectancy), but the higher the variance of the lifetime distribution. The latter aspect is due to the mix of possible pathologies leading to each specific individual classification (and also due to data scarcity). A higher degree of (partially unobservable) heterogeneity in mortality follows, inside each sub-portfolio of special-rate annuities. However, this heterogeneity can be reduced by restricting the range of pathologies that entitle one to a special-rate annuity, then making the relevant sub-portfolio more homogeneous.

It is also worth noting that most of the available mortality statistics refer to specific sets of pathologies (e.g., diabetes), rather than to broad sets of diseases. Splitting a class of special-rate annuities (for example, the impaired life annuities) into pathology-related subclasses can, hence, be an appropriate choice.

To describe the age patterns of mortality in quantitative terms, appropriate life tables can be chosen. However, the use of a mortality law (in particular a "simple" law) significantly eases the implementation of a sensitivity analysis, which can be performed by assigning diverse values to the relevant parameters. In line with the main purpose of life annuities, that is, providing a post-retirement income, adult and old ages have only been addressed in what follows. Then, the Gompertz law has been used to express the age-pattern of mortality. Parameter values have been chosen to represent the features of the curves of deaths, as described in Section 4.1.

In terms of the force of mortality, *µx*, the Gompertz law is as follows:

$$
\mu\_{\mathbf{x}} = B \, \mathbf{c}^{\mathbf{x}}, \quad \text{with } B, \mathbf{c} > \mathbf{0} \tag{1}
$$

Instead of referring to the usual parametrization in (1), we refer to the "informative" parametrization (see, for example, Carriere 1992), that is:

$$
\mu\_{\mathbf{x}} = \frac{1}{D} \exp\left(\frac{\mathbf{x} - M}{D}\right), \quad \text{with } M, D > 0 \tag{2}
$$

where *M* denotes the mode of the Gompertz probability density function and *D* a measure of dispersion. Relations with the usual parameters are as follows:

$$c = \exp\left(\frac{1}{D}\right)\_\text{\,\,\,\,}\tag{3}$$

$$B = \frac{\exp\left(-\frac{M}{D}\right)}{D} \tag{4}$$

Sensitivity analysis can simply be performed by assigning values to the mode parameter *M* and the dispersion parameter *D* (as described in Sections 6.2 and 7.2).

It can be proved that the elements of the corresponding life table {`*x*}, with `<sup>0</sup> = 100,000, are given by the following expression:

$$\ell\_{\rm x} = 100 \beta 000 \times \exp\left(\exp\left(-\frac{M}{D}\right) - \exp\left(\frac{\chi - M}{D}\right)\right) \tag{5}$$

From the life table {`*x*}, all the biometric functions of interest (e.g., *dx*, *qx*, etc.) can immediately be derived.

As already noted, the shape of the lifetime distribution can be driven by choosing specific values *M<sup>k</sup>* and *D<sup>k</sup>* for the parameters of Equation (2). The (baseline) parameter values shown in Table 1 determine the curves of death plotted in Figure 1. We note that the choice of the parameter values is only aimed at performing a sensitivity analysis and does not reflect real statistical data.

**Table 1.** Parameters of the Gompertz law.


#### **5. The Actuarial Model**

After defining the portfolio structures used in the various evaluations, we define the quantities referred to in the deterministic and the stochastic assessments.

*5.1. Portfolio Structures*

We will consider a life annuity portfolio P generally consisting of three sub-portfolios:


Of course, one of the *n<sup>k</sup>* can be set equal to 0. Let *n* denote the size of the portfolio P, that is:

$$n = n\_1 + n\_2 + n\_3$$

Assumptions underlying the actuarial model are as follows:


#### *5.2. Actuarial Values*

Our ultimate object is to analyze the behavior of various quantities defined as functions of *n*1, *n*2, *n*3, in particular: expected value, variance and coefficient of variation (risk index) of the portfolio payouts.

To this purpose, we first recall the basic formulae for a life annuity-immediate, with benefit *b* = 1 paid to an individual age *x* at policy issue and assigned to sub-portfolio SP*k*.

The expected present value (shortly, the actuarial value) of the annual benefits paid to the individual is given (according to the traditional actuarial notation) by:

$$a\_{\mathbf{x}}^{(k)} = \sum\_{h=1}^{\omega - \chi} a\_{h|\lceil h \rceil 1} q\_{\mathbf{x}}^{(k)} \tag{6}$$

where:


For example, with the parameter values given in Table 1 and *x* = 65, Equation (6) yields:

$$\begin{aligned} a\_{65}^{(1)} &= 17.29\\ a\_{65}^{(2)} &= 11.00\\ a\_{65}^{(3)} &= 8.20 \end{aligned}$$

The variance of the present value of the annual benefits is given by:

$$
\sigma\_{\mathbf{x}}^{2(k)} = \sum\_{h=1}^{\omega - \chi} a\_{h\parallel h \mid 1}^2 q\_{\mathbf{x}}^{(k)} - \left(a\_{\mathbf{x}}^{(k)}\right)^2 \tag{7}
$$

For example, with the above data, from Equation (7) we obtain:

$$
\begin{aligned}
\sigma\_{65}^{2(1)} &= 16.858 \\
\sigma\_{65}^{2(2)} &= 26.436 \\
\sigma\_{65}^{2(3)} &= 27.446
\end{aligned}
$$

We denote with E (*k*) (*n<sup>k</sup>* ) and Var(*k*) (*n<sup>k</sup>* ) the expected value and the variance of the benefit payouts of sub-portfolio SP*k*. Assuming a benefit *b* = 1, we obviously have:

$$\mathbb{E}^{(k)}(n\_k) = n\_k a\_x^{(k)}; \ k = 1, 2, 3 \tag{8}$$

and, thanks to the assumption of independence among the individual lifetimes:

$$\text{Var}^{(k)}(n\_k) = n\_k \sigma\_\mathbf{x}^{\mathcal{Z}(k)}; \ k = 1, 2, 3 \tag{9}$$

For a generic portfolio P, consisting of *n* = *n*<sup>1</sup> + *n*<sup>2</sup> + *n*<sup>3</sup> policies, we then find:

$$\mathbb{E}(n\_1, n\_2, n\_3) = \sum\_{k=1}^{3} \mathbb{E}^{(k)}(n\_k) \tag{10}$$

$$\mathbb{Var}(n\_1, n\_2, n\_3) = \sum\_{k=1}^3 \mathbb{Var}^{(k)}(n\_k) \tag{11}$$

#### *5.3. The Risk Index*

The risk index (or coefficient of variation) is a relative risk measure that expresses the variability of a random quantity in terms of standard deviation per unit of expected value. It is frequently adopted in risk theory and risk management to assess the so-called pooling effect, that is, the diversification effect which is achieved by constructing a pool of risks.

For a generic portfolio P, the risk index *ρ* is defined as follows:

$$\rho(n\_1, n\_2, n\_3) = \frac{\sqrt{\text{Var}(n\_1, n\_2, n\_3)}}{\mathbb{E}(n\_1, n\_2, n\_3)}\tag{12}$$

We note that *ρ* is a unit-free risk measure.

#### *5.4. Cash Flows*

Annual cash flows are, of course, random quantities. For the generic sub-portfolio SP*k*, the random cash flow, that is the sub-portfolio payout, at time *t*, *X<sup>k</sup>* (*t*), depends on the number *N<sup>k</sup>* (*t*) of annuitants alive at that time (out of the initial *n<sup>k</sup>* ), and is of course given by:

$$X\_k(t) = b \, N\_k(t);\ \; k = 1, 2, 3 \tag{13}$$

Referring to a generic portfolio P, consisting of three sub-portfolios, the total payout at time *t* is then given by:

$$X(t) = \sum\_{k=1}^{3} X\_k(t) \tag{14}$$

#### **6. Portfolio Risk Profiles: Deterministic Approach**

In this Section, we first assess the impact, in terms of the risk index, of the portfolio structure on the portfolio risk profile. Biometric assumptions are as specified in Section 4.2, with parameter values given in Table 1 (if not otherwise stated).

We then assess the impact, again in terms of the risk index, of diverse biometric assumptions.

#### *6.1. Impact of the Portfolio Structure*

#### 6.1.1. Cases 1.1

We analyze the impact of the size of the sub-portfolio SP2 of enhanced annuities. Then:

$$\begin{aligned} n\_1 &= 10 \text{\AA} \\ n\_2 &= 100 \text{\AA} \\ n\_3 &= 0 \end{aligned}$$

Results are shown in Table 2.


**Table 2.** Cases 1.1—Impact of the portfolio structure on the risk index.

#### 6.1.2. Cases 1.2

We analyze the impact of the size of the sub-portfolio SP3 of impaired annuities. Then:

$$\begin{aligned} n\_1 &= 10 \text{\AA} \\ n\_2 &= 0 \\ n\_3 &= 100 \text{\AA} 200 \dots \text{\AA} \end{aligned}$$

Results are shown in Table 3.



#### 6.1.3. Cases 1.3

We assume that both enhanced annuities and impaired annuities are sold (together with standard annuities), and analyze the joint impact by assuming that *n*<sup>3</sup> = *n*2/2. Then:

$$\begin{aligned} n\_1 &= 10 \text{\AA} \\ n\_2 &= 500, 600, \dots, 1000 \\ n\_3 &= 250, 300, \dots, 500 \end{aligned}$$

Results are shown in Table 4.


**Table 4.** Cases 1.3—Impact of the portfolio structure on the risk index.

#### 6.1.4. Cases 1.4

The launch of special-rate annuities might negatively impact on the sale of standard annuities (the so called "cannibalization effect"). To analyze this aspect in terms of portfolio risk profile, we assume that one half of the enhanced annuity sales (sub-portfolio SP2) are "subtracted" from the standard annuity business (sub-portfolio SP1). Then, we consider portfolios with the following sub-portfolio sizes:

$$\begin{aligned} n\_1 &= 10,000 - \frac{n\_2}{2} \\ n\_2 &= 500,600 \dots, 1000 \\ n\_3 &= \frac{n\_2}{2} = 250,300 \dots, 500 \end{aligned}$$

Furthermore, it is reasonable to assume that, in case of a cannibalization effect, the mortality in the standard annuity sub-portfolio improves. To represent this aspect, we assume *M*<sup>1</sup> = 91 (instead of *M*<sup>1</sup> = 90). Results are shown in Table 5.


**Table 5.** Cases 1.4—Impact of the portfolio structure on the risk index.

#### 6.1.5. Some Comments

When considering a given set of cases, the size of subportfolios and structure of the total portfolio change, and this of course impacts both the numerator and denominator of the risk index. Hence, the analysis of the risk index values in the various portfolio structures provides interesting information. We note that, in all the sets of cases we have considered, the range of values assumed by the risk index is very narrow. From a mathematical perspective, this is the straight consequence of a higher variability in terms of standard deviation (the numerator of fraction (12)) offset, to a large extent, by a higher expected value (the denominator), and, in practice, a higher volume of premiums. Therefore, the almost constant value of the risk index witnesses this offset. A wider range of values (anyway very limited) can be noted as the effect of the number of impaired annuities: see, for example, the set of cases 1.2, where the increase in the risk index is equal to 2.6%, compared to the set 1.1 where the increase is smaller than 1%.

The presence of special-rate annuities impacts on the standard annuity portfolios, in particular, in terms of the lifetime distribution of standard annuitants. This has been considered in Section 6.1.4 by increasing the modal age at death from 90 to 91, then increasing the expected value of benefits paid by standard annuities. It is worth noting that a (reasonable) impact on standard annuity premiums should follow, and this might, in turn, impact the demand of standard annuities. Further interesting results, regarding the

variability of the annual payouts, can be achieved via stochastic analysis and are presented in Section 7.

#### *6.2. Impact of Lifetime Distributions*

Given the uncertainty in biometric assumptions, a sensitivity analysis is appropriate. While keeping unchanged the parameters *M<sup>k</sup>* (representing the modal age at death), we propose diverse assumptions regarding the dispersion of the lifetime distributions, which might more heavily impact on the portfolio risk profile. Hence, various values of the parameters *D<sup>k</sup>* are considered.

Figure 2 shows the graphs of the lifetime distribution for enhanced annuities, corresponding to different values of dispersion (parameter *D*) while keeping the same modal value (parameter *M*).

**Figure 2.** Three assumptions on lifetime dispersion for enhanced annuities.

#### 6.2.1. Cases 2.1

We consider a portfolio only consisting of standard annuities and enhanced annuities, with given sub-portfolio sizes. Hence:

$$\begin{aligned} n\_1 &= 10 \text{J000} \\ n\_2 &= 1000 \\ n\_3 &= 0 \end{aligned}$$

We analyze the impact of diverse assumptions on the dispersion of lifetimes in sub-portfolio SP2. Then:

$$D\_2 = 4, 5, \dots, 13$$

(while keeping *D*<sup>1</sup> = 5). Results are shown in Table 6.


**Table 6.** Cases 2.1—Impact of the lifetime distribution on the risk index.

#### 6.2.2. Cases 2.2

We consider a portfolio only consisting of standard annuities and impaired annuities, with given sub-portfolio sizes. Hence:

$$\begin{aligned} n\_1 &= 10 \, 000 \\ n\_2 &= 0 \\ n\_3 &= 1000 \end{aligned}$$

We analyze the impact of diverse assumptions on the dispersion of lifetimes in sub-portfolio SP3. Then:

$$D\_3 = 11, 12, \dots, 15$$

Results are shown in Table 7.

**Table 7.** Cases 2.2—Impact of the lifetime distribution on the risk index.


6.2.3. Cases 2.3

We consider a portfolio consisting of standard annuities, enhanced annuities and impaired annuities, with given sub-portfolio sizes. Hence:

> *n*<sup>1</sup> = 10,000 *n*<sup>2</sup> = 1000 *n*<sup>3</sup> = 500

We analyze the joint impact of diverse assumptions on the dispersion of lifetimes in both sub-portfolios SP2 and SP3. To this purpose, we assume:

$$D\_2 = D\_3 = 4, 5, \dots, 13$$

We note that lower dispersions can be achieved by restricting the range of pathologies, which entitle the purchase of enhanced annuities and impaired annuities. Results are shown in Table 8.


**Table 8.** Cases 2.3—Impact of the lifetime distributions on the risk index.

#### 6.2.4. Some Comments

Although dispersion in lifetime distributions does affect the risk profile of the annuity portfolio, the sensitivity analysis we have performed witnesses a rather limited impact on the risk index. We note that, of course, the broadest range of risk index values can be found when a portfolio consisting of standard annuities, enhanced annuities and impaired annuities is addressed, and for both the types of special-rate annuities higher values for the dispersion parameter are considered.

#### **7. Portfolio Risk Profiles: Stochastic Approach**

Deterministic assessments performed in Section 6 only provide values of specific markers, notably the risk index. To obtain better insights into the risk profile of a portfolio, stochastic assessments are required. To this purpose, stochastic (Monte Carlo) simulation procedures are commonly adopted.

As we focus on the biometric features of the various portfolios, simulation of the numbers of survivors, that is *N<sup>k</sup>* (*t*), *k* = 1, 2, 3 and *t* = 1, 2, . . . , is only needed. Then, via Equations (13) and (14), the simulated outcomes of the payouts and, finally, the relevant (empirical) distributions are obtained. Consistent with the approach adopted in Section 6, we assume that all the assessments are performed at time *t* = 0, and hence our information is given by the initial sizes of the sub-portfolios, that is, *n*1, *n*<sup>2</sup> and *n*3.

The following Sections 7.1 and 7.2 are organized similarly to Sections 6.1 and 6.2, respectively, but with a reduction in the number of cases analyzed.

Besides "descriptive" results in terms of (empirical) distributions of the annual payouts, the stochastic approach can also yield "operational" results: an example is provided in Section 7.3, where amounts of assets are calculated, which are needed to meet the annual payouts with an assigned probability.

#### *7.1. Impact of the Portfolio Structure*

As already noted, we follow the organization in the cases adopted in Section 6.1, although reducing the number of alternatives.

7.1.1. Cases 1.1

We analyze the impact of the size of the sub-portfolio SP2 of enhanced annuities. Then:

$$\begin{aligned} n\_1 &= 10,000 \\ n\_2 &= 100,500,1000 \\ n\_3 &= 0 \end{aligned}$$

Empirical distributions at times 5 and 10 of the portfolio payout, that is, empirical distributions of *X*(5) and *X*(10), are sketched in Figures 3 and 4, where the three portfolios are denoted by P01, P02 and P03, respectively.

**Figure 3.** Impact of the number of enhanced annuities: empirical distributions of the annual benefit payout at time *t* = 5.

**Figure 4.** Impact of the number of enhanced annuities: empirical distributions of the annual benefit payout at time *t* = 10.

7.1.2. Cases 1.4

To assess the impact of a possible cannibalization effect, we consider three portfolios with the following sub-portfolio sizes:

$$\begin{aligned} n\_1 &= 10,000 - \frac{n\_2}{2} \\ n\_2 &= 500,800,1000 \\ n\_3 &= \frac{n\_2}{2} = 250,400,500 \end{aligned}$$

As previously noted, to represent an improvement in mortality in the standard annuity sub-portfolio, we assume *M*<sup>1</sup> = 91 (instead of *M*<sup>1</sup> = 90). Empirical distributions of the portfolio payout *X*(5) and *X*(10) are sketched in Figures 5 and 6, respectively, where the three portolios are denoted by P01, P02 and P03.

**Figure 5.** Impact of cannibalization effect: empirical distributions of the annual benefit payout at time *t* = 5.

**Figure 6.** Impact of cannibalization effect: empirical distributions of the annual benefit payout at time *t* = 10.

#### 7.1.3. Some Comments

Results are self-explanatory, and in line with the findings in the deterministic setting: the larger the (initial) number of enhanced annuities in Cases 1.1, the higher the dispersion in the annual payouts. The same effect is, of course, witnessed by the distributions of payouts in Cases 1.4, where presence of both enhanced annuities and impaired annuities is assumed.

#### *7.2. Impact of the Lifetime Distribution*

To assess the impact of uncertainty in biometric assumptions, we only analyze the Cases 2.1 considered in Section 6.2.

#### 7.2.1. Cases 2.1

We consider a portfolio only consisting of standard annuities and enhanced annuities, with given sub-portfolio sizes. Hence:

$$\begin{aligned} n\_1 &= 10 \text{\AA} \\ n\_2 &= 1000 \\ n\_3 &= 0 \end{aligned}$$

We analyze the impact of diverse assumptions on the dispersion of lifetimes. Then:

$$D\_2 = 4, \ 6, \ 8, \ 10, \ 12$$

(while keeping *D*<sup>1</sup> = 5). Empirical distributions of the portfolio payout *X*(5) and *X*(10) are sketched in Figures 7 and 8, where the portfolios are respectively denoted by P01, . . . , P05.

**Figure 7.** Impact of different dispersion parameters: empirical distributions of the total payout at time *t* = 5.

**Figure 8.** Impact of different dispersion parameters: empirical distributions of the total payout at time *t* = 10.

#### 7.2.2. Some Comments

These results are self-evident: a larger dispersion in the lifetime distribution of annuitants with enhanced annuity implies a larger dispersion in the total portfolio payout (compare, in particular, distributions in portfolio P05 and P01. Moreover, for all the portfolios, dispersions increase with time (compare the distributions in Figure 7 to the ones in Figure 8).

#### *7.3. Meeting the Annual Payouts*

Appropriate resources must be assigned to the portfolio in order to meet the annual payouts with a high probability. Diverse criteria can be adopted to quantify the above resources which, whatever the criterion adopted, will partly be provided by single premiums cashed at policies issued (via decumulation of the portfolio reserve) and partly by shareholders' capital allocated to the portfolio. In what follows, we focus on the annual total amount of resources needed, disregarding the funding source.

#### 7.3.1. The Percentile Principle

Referring to a generic portfolio and the relevant cash flows, we recall that *X*1(*t*), *X*2(*t*), *X*3(*t*) denote the random payouts at time *t*, related to standard annuities, enhanced annuities and impaired annuities, respectively (see Section 5.4), and:

$$X(t) = X\_1(t) + X\_2(t) + X\_3(t)\tag{15}$$

denotes the portfolio total payout at time *t*.

We adopt the percentile principle. Hence, we have to find, for *t* = 1, 2, . . . , the amount *A*(*t*) such that:

$$\Pr[X(t) > A(t)] = \varepsilon \tag{16}$$

where *ε* denotes an assigned (small) probability; hence, 1 − *ε* can be interpreted as the "adequacy" level.

A more detailed analysis could be performed by separately addressing the risk profile of each sub-portfolio, thus calculating, for *k* = 1, 2, 3 and *t* = 1, 2, . . . the quantities *A<sup>k</sup>* (*t*) such that:

$$\Pr[X\_k(t) > A\_k(t)] = \varepsilon\_k \tag{17}$$

However, we only focus on the overall requirement (16), which clearly takes into account the pooling effect.

#### 7.3.2. Numerical Results

We consider four portfolios with the structures defined in Table 9.

#### **Table 9.** Portfolio structures.


This way, we can analyze the risk profile of a "traditional" portfolio only consisting of standard annuities (P01), a portfolio including standard annuities and enhanced annuities (P02) or standard annuities and impaired annuities (P03) and finally a portfolio including both types of special-rate annuities (P04).

Values of the parameters *M<sup>k</sup>* and *D<sup>k</sup>* for standard life annuities (*k* = 1), enhanced annuities (*k* = 2) and impaired annuities (*k* = 3), respectively, are as specified in Table 1.

Asset requirements (at a given time *<sup>t</sup>*), in terms of ratios Assets *EV* , where *EV* denotes the expected value of the total payout, are plotted in Figures 9–11 against the probability 1 − *ε*.

**Figure 9.** Impact of the portfolio structure on assets requirements at time *t* = 1.

**Figure 10.** Impact of the portfolio structure on assets requirements at time *t* = 5.

**Figure 11.** Impact of the portfolio structure on assets requirements at time *t* = 10.

#### 7.3.3. Some Comments

Results in terms of assets requirements are also encouraging. We note that the range of values, expressed by the ratio between assets required and expected values of total payout, corresponding to the various portfolio structures are very limited, whatever the adequacy level chosen. A higher sensitivity with respect to the adequacy level 1 − *ε* can be observed in particular for *t* = 10, because of a dispersion of the numbers of survivors and hence of the total payouts, which increases with time.

We note that the portfolio structure of course evolves over time and, notably, the share of standard annuities progressively increases because of longer lifetimes of standard annuitants. Increasing shares of standard life annuities (characterized, according to a reasonable assumption, by a variance lower than that of special-rate annuities) imply a lower riskiness in the (total) annuity portfolio. In terms of the ratio Assets/EV, this effect appears in particular when comparing the requirements (according to each probability 1 − *ε*) represented in Figures 9–11.

#### **8. Concluding Remarks**

By offering special-rate life annuities, on the one hand, a higher premium income can be expected, while on the other hand, a higher variability of the total portfolio payout will follow because of both the larger size and the specific higher variability of payouts related to special-rate annuities.

The analysis, in quantitative terms, of the "balance" between the two aspects (that is, higher risk and higher premium income) has been the aim of this research. A number of numerical evaluations have been performed by adopting both a deterministic approach and a stochastic one as well. Diverse hypotheses on lifetime distributions have been assumed, and various portfolio sizes and structures (in terms of numbers of standard, enhanced and impaired annuities) have been considered.

It is worth noting that, whatever the choice of the parameter values for the Gompertz law, the mortality model is deterministic, in the sense that no uncertainty in future mortality trend is accounted for. In this context, it is reasonable to assume independence among the individual lifetimes. An interesting extension of our simple model should allow for uncertainty in mortality trends via an appropriate stochastic mortality model. Then, correlation among lifetimes would follow, so that independence assumption would be replaced by conditional independence. Moreover, diverse trends and diverse degrees of uncertainty could be considered for the various special-rate annuities by taking into account possible improvements in medical treatments, surgery, etc.

Special attention should be placed on the "cannibalization" effect. First, a reduction in the size of the standard annuity sub-portfolio may occur, provided that some applicants can be eligible to special-rate annuities. A lower average mortality in the standard sub-portfolio then follows (as noted in Section 6.1.5), leading to a (reasonable) increase in standard annuity premiums. This might in turn impact the demand of standard life annuities (even beyond the reduction mentioned above).

A further aspect that is interesting to investigate is the impact on the insurer's liabilities (and, notably, on the portfolio risk profile) of an incorrect allocation of individual risks to the various rating classes. Because of the presence of unobservable risk factors, a misspecification of the rating class is always possible. This would result in an unfair annuity rate applied to some individuals and then in an increased (or reduced) probability of loss for the insurer. In this regard, the research task should concern, in particular, the modeling of the incorrect specification of the rating class.

The results we have obtained of course depend on assumptions (notably, regarding both the portfolio structure, the mortality model and the relevant parameters). Nevertheless, the broad range of assumptions regarding both the portfolio structure and the lifetime distributions has allowed us to perform an effective sensitivity analysis, whose interesting achievements witness the possibility of extending the life annuity business without taking huge amounts of risk. Hence, the creation of values for customers (and an increase in the insurer's market share) can be pursued without a significant worsening of the company's risk profile.

**Author Contributions:** Conceptualization, E.P. and D.Y.T.; methodology, E.P. and D.Y.T.; software, E.P. and D.Y.T.; validation, E.P. and D.Y.T.; formal analysis, E.P. and D.Y.T.; investigation, E.P. and D.Y.T.; writing—original draft preparation, E.P. and D.Y.T.; writing—review and editing, E.P. and D.Y.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors thank the anonymous referees for their helpful comments and suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Note**

<sup>1</sup> Stage-Gate® is a registered trademark of Stage-Gate Inc. (www.stage-gate.com (accessed on 10 January 2022)).

#### **References**

Ainslie, Ross. 2000. *Annuity and Insurance Products for Impaired Lives*. Working Paper. London: Staple Inn Actuarial Society.

Beard, Robert E. 1959. Note on some mathematical mortality models. In *CIBA Foundation Colloquia on Ageing*. Edited by Gordon Ethelbert Ward Wolstenholme and Maeve O. Connor. Boston: John Wiley & Sons, Ltd., vol. 5, pp. 302–11.

Carriere, Jacques F. 1992. Parametric models for life tables. *Transaction of Society of Actuaries* 44: 77–99.

Denuit, Michel, and Esther Frostig. 2006. Heterogeneity and the need for capital in the individual model. *Scandinavian Actuarial Journal* 2006: 42–66. [CrossRef]

Denuit, Michel, and Esther Frostig. 2007. Association and heterogeneity of insured lifetimes in the Lee–Carter framework. *Scandinavian Actuarial Journal* 2007: 1–19. [CrossRef]

Drinkwater, Matthew, Joseph E. Montminy, Eric T. Sondergeld, Christopher G. Raham, and Chad R. Runchey. 2006. Substandard Annuities. Technical Report. LIMRA International Inc. and the Society of Actuaries, in Collaboration with Ernst & Young LLP. Available online: https://www.soa.org/Files/Research/007289-Substandard-annuities-full-rpt-REV-8-21.pdf (accessed on 10 January 2022).


Gatzert, Nadine, Gudrun Schmitt-Hoermann, and Hato Schmeiser. 2012. Optimal risk classification with an application to substandard annuities. *North American Actuarial Journal* 16: 462–86. [CrossRef]

Gracie, Stewart, and Stephen Makin. 2006. The Price to Pay for Enhanced Annuities. Healthcare Conference 2006. Available online: https://www.actuaries.org.uk/system/files/documents/pdf/Gracie.pdf (accessed on 10 January 2022).

Haberman, Steven, and Annamaria Olivieri. 2014. Risk Classification/Life. In *Wiley StatsRef: Statistics Reference Online*. Hoboken: Wiley. Hoermann, Gudrun, and Jochen Russ. 2008. Enhanced annuities and the impact of individual underwriting on an insurer's profit

situation. *Insurance: Mathematics & Economics* 43: 150–57.

James, Mick. 2016. Enhanced annuities: Caring for at-retirement needs. *Reinsurance News* 84: 24–27.


Pitacco, Ermanno. 2019. Heterogeneity in mortality: A survey with an actuarial focus. *European Actuarial Journal* 9: 3–30. [CrossRef]

Pitacco, Ermanno. 2020. *ERM and QRM in Life Insurance. An Actuarial Primer*. Springer Actuarial. Berlin: Springer.

Pitacco, Ermanno. 2021. *Life Annuities. From Basic Actuarial Models to Innovative Designs*. London: Risk Books.


## *Article* **The Copula Derived from the SAHARA Utility Function**

**Jaap Spreeuw**

Faculty of Actuarial Science and Insurance, Bayes Business School (Formerly Cass), University of London, 106 Bunhill Row, London EC1Y 8TZ, UK; j.spreeuw@city.ac.uk

**Abstract:** A new Archimedean copula family is presented that was derived from the SAHARA utility function introduced in the economic literature in 2011. Its properties are discussed, and its flexibility and versatility are demonstrated. It is left tail decreasing or right tail increasing, but unlike mainstream Archimedean families, not necessarily stochastically increasing at the same time. It is shown that the family fits very well to a dataset of previously studied coupled lives in the literature.

**Keywords:** copula; Archimedean generator; dependence; coupled lives

#### **1. Introduction**

Archimedean copulas, which in two dimensions are of the form *Cψ*(*u*1, *u*2) = *ψ ψ* −1 (*u*1) + *ψ* −1 (*u*2) where *ψ* is the one-dimensional generator, have become a popular mode of modelling dependence in both finance and insurance. Several ways of constructing copula families are given in Chapter 3 of Joe (2015). The interpretation of the generator as the Williamson transform of a radial random variable has given rise to new Archimedean families; see McNeil and Nešlehová (2009, 2010). Archimedean copulas are a flexible class due to the ease with which new Archimedean copulas with an enriched parameter space can be constructed from existing ones using transformations. For the bivariate case, five of such transformations, namely, left composition, right composition, scaling, exponentiation, and the linear combination of the (inverse) generator, were introduced in the literature by Genest et al. (1998). They were reviewed by Michiels and De Schepper (2012), and in more detail in Michiels and De Schepper (2009), with the focus on the so-called *λ* function (which is the ratio of the inverse generator to its derivative). For the latter, see also Michiels et al. (2011).

In the literature, several generalized families were constructed that contain the Archimedean class as a special case, e.g., the Archimax family in Capéraà et al. (2000) (which includes extreme value copulas as another special case). The background risk model where one random variable ("systemic risk") acts multiplicatively on a series of other random variables ("idiosyncratic risks") is the basis of generalization of Archimedean copulas in several ways, as demonstrated in Côté and Genest (2019) and Marri and Moutanabbir (2022).

Commonly used families typically feature a generator being completely monotonic and thereby the Laplace transform of a mixing random variable. Common examples include Clayton and Frank (as far as dependence is positive), Gumbel–Hougaard, and Joe. This subfamily of Archimedean copulas, also known as shared frailty models, has the advantage of being valid in any dimension. Recent applications regarding shared frailty models involve the aforementioned background risk model; see (Albrecher et al. 2011; Furman et al. 2021; Sarabia et al. 2018).

The Gumbel–Hougaard copula is a good fit to the well-known dataset of loss vs. Allocated Loss Adjustment Expenses (ALAE), which is the object of statistical inference in several publications, starting with Genest et al. (1998). For extensive analysis, consult Joe (2015). Probably there are several other case studies of dependence in insurance where the use of shared frailty models is appropriate.

**Citation:** Spreeuw, Jaap. 2022. The Copula Derived from the SAHARA Utility Function. *Risks* 10: 133. https://doi.org/10.3390/ risks10070133

Academic Editors: Ermanno Pitacco and Annamaria Olivieri

Received: 28 February 2022 Accepted: 10 June 2022 Published: 28 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Shared frailty dependency models have the property of being conditionally increasing (CI) and multivariate totally positive of order 2 (MTP2), as shown in Müller and Scarsini (2005). In the specific case of two dimensions, this is known as TP2, implying stochastically increasing (SI), which in turn implies both left tail decreasing (LTD) and right tail increasing (RTI). If the two involved marginal random variables are remaining lifetimes, LTD/RTI implies that the hazard rate of the one upon hazard (e.g., death or default) of the other goes up. The stronger condition SI implies that the hazard rate of the one given the hazard of the other (e.g., death or default) at time *t* decreases as *t* goes up. This was noted in Spreeuw (2006). The underlying assumption of SI may work out well in reliability theory. Consider two printers available to office staff. When one fails, the other one, pending the repair of the first, is used more intensively and thereby exposed to greater strain. It is sensible to assume that the longer the first printer is out of order, the higher the failure rate of the remaining machine would be.

For coupled lives, however, it is not so clear-cut. On the one hand, the event of the death of one life usually triggers an elevated mortality of the surviving life, so the assumption of LTD/RTI seems sensible. In addition, such lives are exposed to common risks due to permanently living together and due to a similar background ("birds of a feather flock together", the so called long-term dependence according to Hougaard (2000)). On the other hand, however, there is also the phenomenon of the event of the one life dying leading to the mortality of the remaining life temporarily going up, which is the so-called broken-heart syndrome (short-term dependence, also attributed to Hougaard (2000)). Similar nonstandard features may apply to other cases in insurance (and finance as well). In short, there is a case for constructing copula families allowing for flexibility in terms of type of dependence, such as LTD/RTI but not necessarily SI.

In this paper, we introduce a new Archimedean copula family that is based on a link between Archimedean generators and utility functions; see Spreeuw (2010) for more details. Unlike mainstream copulas, this family has the property of being LTD/RTI, but not necessarily SI, the latter being clearly indicated by the sign of one of the parameters.

The outline of the paper is as follows. Section 2 gives the basic definitions of Archimedean copula, and the dependence concepts of LTD/RTI and SI. Section 3 introduces the new family and analyzes its basic properties. Section 4 fits the new Archimedean family to the section of censored remaining lifetime data of coupled lives, as in Luciano et al. (2008). Section 5 sets out a conclusion.

#### **2. Basic Definitions**

Define *ψ* as the generator of a 2-dimensional Archimedean copula, being strictly continuous, strictly decreasing, convex, with *ψ*(0) = 1 and lim*x*→<sup>∞</sup> *ψ*(*x*) = 0. The copula itself is then specified as

$$\mathbb{C}\_{\psi}(\boldsymbol{u}\_1, \boldsymbol{u}\_2) = \psi\Big(\psi^{-1}(\boldsymbol{u}\_1) + \psi^{-1}(\boldsymbol{u}\_2)\Big),\tag{1}$$

where *u*1, *u*<sup>2</sup> each take values between 0 and 1.

Next are definitions of the tail concepts of left tail decreasing (LTD), right tail increasing (RTI) and stochastically increasing (SI), based on two random variables *X* and *Y*, and their copula *C*. They can all be found in Chapter 5 of Nelsen (2006).

**Definition 1.** *Y is LTD in X (notation LTD*(*Y*|*X* )*)* ⇔ Pr[*Y* ≤ *y*|*X* ≤ *x* ] *is nonincreasing in x for all y. For an exchangeable copula C (i.e., C*(*u*, *v*) = *C*(*v*, *u*) *for* 0 ≤ *u*, *v* ≤ 1*) of random variables X and Y, LTD*(*Y*|*X* ) *and LTD*(*X*|*Y*) *are equivalent.*

**Definition 2.** *Y is RTI in X (notation RTI*(*Y*|*X* )*)* ⇔ Pr[*Y* > *y*|*X* > *x* ] *is nondecreasing in x for all y. For an exchangeable copula C (i.e., C*(*u*, *v*) = *C*(*v*, *u*) *for* 0 ≤ *u*, *v* ≤ 1*) of random variables X and Y, RTI*(*Y*|*X* ) *and RTI*(*X*|*Y*) *are equivalent.*

**Definition 3.** *Y is SI in X (notation SI*(*Y*|*X* )*)* ⇔ Pr[*Y* ≤ *y*|*X* = *x* ] *is nonincreasing in x for all y*.

The following propositions are from Avérous and Dortet-Bernadet (2004). The second one was originally shown in Capéraà and Genest (1993).

**Proposition 1.** *If C is Archimedean with generator ψ, LTD*(*Y*|*X* ) *or LTD*(*X*|*Y*) *if and only if ψ is logconvex. Likewise, if C*ˆ*, the rotated copula (also known as survival copula) of C (so <sup>C</sup>*ˆ(*u*, *<sup>v</sup>*) <sup>=</sup> *<sup>u</sup>* <sup>+</sup> *<sup>v</sup>* <sup>−</sup> <sup>1</sup> <sup>+</sup> *<sup>C</sup>*(<sup>1</sup> <sup>−</sup> *<sup>u</sup>*, 1 <sup>−</sup> *<sup>v</sup>*) *for* <sup>0</sup> <sup>≤</sup> *<sup>u</sup>*, *<sup>v</sup>* <sup>≤</sup> <sup>1</sup>*) is Archimedean with generator <sup>ψ</sup>*<sup>b</sup> *is RTI*(*Y*|*X* ) *or RTI*(*X*|*Y*) *if and only if ψ*b *is logconvex.*

**Proposition 2.** *If C is Archimedean with differentiable generator ψ, SI*(*Y*|*X* ) *or SI*(*X*|*Y*) *if and only if* −*ψ* 0 *is logconvex.*

#### **3. SAHARA Family**

The SAHARA copula family is derived from the Symmetric Asymptotic Hyperbolic Absolute Risk Aversion (SAHARA) utility function introduced in Chen et al. (2011). This utility function is specified below.

$$\varrho\_{\theta,\varepsilon}(s)$$

$$
\mathcal{J} = \begin{cases}
\cdot \begin{pmatrix} s - \varepsilon + \left(1 + 1/\theta\right)\sqrt{\delta^2 + \left(s - \varepsilon\right)^2} \\
\frac{1}{2}\ln\left(\left(s - \varepsilon\right) + \sqrt{\delta^2 + \left(s - \varepsilon\right)^2}\right) \\
\end{pmatrix} & \theta \in \left(-\infty, -1\right) \cup \left(0, \infty\right), \delta > 0, \varepsilon \in \mathbb{R}. \tag{2}
$$

As shown in Spreeuw (2010), a strict Archimedean generator can be obtained from a utility function *ϕ* if *ϕ*(∞) = lim*s*−→<sup>∞</sup> *ϕ*(*s*) < ∞. For SAHARA, this is the case when *θ* > 0. Then, applying the formula *ψθ*,*e*(*s*) = {*ϕθ*,*e*(∞) − *ϕθ*,*e*(*s*)}{*ϕθ*,*e*(∞) − *ϕθ*,*e*(0)} −1 leads to the Archimedean generator

$$\psi\_{\theta,\varepsilon}(s) = \left(\frac{s-\varepsilon+\sqrt{\delta^2+\left(s-\varepsilon\right)^2}}{-\varepsilon+\sqrt{\delta^2+\varepsilon^2}}\right)^{-(1+1/\theta)} \left(\frac{s-\varepsilon+(1+1/\theta)\sqrt{\delta^2+\left(s-\varepsilon\right)^2}}{-\varepsilon+(1+1/\theta)\sqrt{\delta^2+\varepsilon^2}}\right). \tag{3}$$

**Remark 1.** *This approach of obtaining an Archimedean generator from a utility function is not to be confused with the method of obtaining the* inverse *of an Archimedean generator from a utility function. For the latter, consult Spreeuw (2014).*

The SAHARA utility function was inspired by nonmonotone risk aversion coefficient

$$AR\_{\theta}(s) = \frac{1 + 1/\theta}{\sqrt{(s - \epsilon)^2 + \delta^2}}\sqrt{$$

which, unlike common utility functions, is not monotone in its argument. It is rather increasing for *s* < *e* attaining a maximum for *s* = *e*, and decreasing for *s* > *e*. The SAHARA utility function found applications in both finance and insurance (Bernard et al. 2021; Bernard and Kwak 2016; Brachetta and Schmidli 2020; Chen et al. 2021; Chen and Vellekoop 2017; Li and Ma 2018; Schumacher 2018). As shown in Spreeuw (2010), a risk aversion monotone decreasing (increasing) on the positive real line in general implies that the corresponding Archimedean copula is stochastic increasing (stochastic decreasing) in two dimensions. The former property applies to the vast majority of commonly applied

copula families (including all those whose generator is a completely monotonic function). For *e* > 0 the copula is neither stochastic increasing nor stochastic decreasing. To the best of our knowledge, there are hardly any copula families that share this property.

For *e* < 0, the condition imposed on *δ* can be somewhat relaxed to *δ* ≥ 0, *δ* = 0 corresponding to the Clayton copula with parameter *θ*. Due to scaling, for some realvalued *δ* ∗ > 0 and *e* ∗ , (*δ*, *e*) = (*δ* ∗ , *e* ∗ ) gives exactly the same Archimedean copula as (*δ*, *e*) = 1, *e* ∗ *δ* ∗ . To see this, we divide in (3) both numerator and denominator by *δ* > 0. This gives:

$$\psi\_{\theta,\varepsilon}(s) = \left(\frac{\frac{s}{\delta} - \frac{\epsilon}{\delta} + \sqrt{1 + \left(\frac{s}{\delta} - \frac{\epsilon}{\delta}\right)^2}}{-\frac{\epsilon}{\delta} + \sqrt{1 + \left(\frac{\epsilon}{\delta}\right)^2}}\right)^{-(1+1/\theta)} \left(\frac{\frac{s}{\delta} - \frac{\epsilon}{\delta} + (1+1/\theta)\sqrt{1 + \left(\frac{s}{\delta} - \frac{\epsilon}{\delta}\right)^2}}{-\frac{\epsilon}{\delta} + (1+1/\theta)\sqrt{1 + \left(\frac{\epsilon}{\delta}\right)^2}}\right).$$

Now, *e* being real-valued implies that *e*/*δ* is real-valued as well. In addition, it is wellknown that a generator is only defined up to a multiple constant. In other words, for *β* > 0, *ψθ*,*e*(*s*) and *ψθ*,*e*(*β s*) generate the same Archimedean copula. *δ* <sup>∗</sup> −→<sup>+</sup> <sup>0</sup> is equivalent to *e* ∗ *δ* ∗  −→ ∞. Hence, without loss of generality, we take *δ* = 1 from now on bearing in mind that, for *e* −→ −∞, the Clayton copula with parameter *θ* is obtained as a limiting case. For *e* −→ ∞, the Clayton copula is again obtained as a limiting case, although now with parameter <sup>−</sup> *<sup>θ</sup>* 2*θ*+1 .

It can be numerically shown that this copula (a) for a fixed *θ* decreases in concordance with increasing *e*; (b) for negative fixed *e*, it increases in concordance with increasing *θ*; and (c) for a positive fixed *e*, it first decreases and then increases in concordance in terms of *θ*. For any finite *e*, *θ* ↓ 0 and *θ* −→ ∞ lead to independence and comonotonicity, respectively.

According to Theorem 4.3 of Joe (1997), p. 91, for Archimedean copulas with a strict generator, the population version of Kendall's tau can be written as:

$$\tau = 1 - 4 \int\_{s=0}^{\infty} s \left( \frac{d}{ds} \psi\_{\theta, \varepsilon}(s) \right)^2 ds.$$

Using software system *Wolfram Mathematica* (see Wolfram Research Inc. (2017)) this gives the expression:

$$\tau = 1 - \frac{(2\theta + 1)\left( (\theta + 2)(3\theta + 2) + 4(\theta + 2)\epsilon^4 + 12(\theta + 1)\epsilon^2 + 2(\theta + 4)\epsilon\sqrt{\epsilon^2 + 1} + 4(\theta + 2)\epsilon^2\sqrt{\epsilon^2 + 1} \right)}{(\theta + 2)(3\theta + 2)\left(\theta + \epsilon\left(\sqrt{\epsilon^2 + 1} + \epsilon\right) + 1\right)^2} \times$$

for *e* = 0 considerably reducing to *τ* = { *θ*/(*θ* + 1)} 2 . Taking the limit for *e* → −∞ , keeping *θ* constant, gives *τ* → *θ*/(*θ* + 2). This is the well-known formula of Kendall's tau for the Clayton copula, and therefore not very surprising. Taking the limit *e* → ∞, keeping *θ* constant, gives *τ* → −*θ*/(3*θ* + 2). This implies that, for increasing *θ*, the range of values taken by *τ* increases. So, the lowest possible value of *τ* for this family is lim*θ*−→<sup>∞</sup> −*θ*/(3*θ* + 2) = −1/3. For fixed nonpositive *e*, *τ* is monotone increasing in *θ* from 0 to 1. For fixed and finite positive *e*, *τ* as a function of *θ* first decreases until a certain negative minimum that is greater than −1/3, and increases afterwards. The greater the value of *e* is, the greater the value of *θ* for which the minimum is reached and the smaller the minimal value. The difference between Kendall's tau for *e* = 0 and *e* → −∞ is

$$\frac{\theta}{\theta+2} - \left(\frac{\theta}{\theta+1}\right)^2 = \frac{\theta}{\left(\theta+2\right)\left(\theta+1\right)^2}.$$

which is zero for *θ* = 0, increasing until *θ* = √ 5 − 1 .<sup>2</sup> <sup>≈</sup> 0.618 (which for *<sup>e</sup>* <sup>=</sup> <sup>0</sup> and *e* → −∞ gives values of Kendall's tau of 0.236 and 0.146, respectively), then decreasing for increasing *θ* and ultimately vanishing. There is no upper tail dependence, while the lower tail dependence coefficient is 2<sup>−</sup> 1/*<sup>θ</sup>* .

Another interesting feature of this family concerns the conditional survival copula, that is, if the copula of the conditional joint survival function *H*(**x**|**y** ) = Pr[*X*<sup>1</sup> > *x*1, *X*<sup>2</sup> > *x*<sup>2</sup> |*X*<sup>1</sup> > *y*1, *X*<sup>2</sup> > *y*<sup>2</sup> ]. If the joint survival function *H*(**x**) has an Archimedean copula with generator *ψ*(*s*),*s* ≥ 0, the conditional joint survival function *H*(**x**|**y** ) also has an Archimedean copula with generator

$$
\psi\_y(s) = \psi(s+t) / \psi(t) \,. \tag{4}
$$

where *t* = *ψ* −1 *H*(**y**) . (Usually the conditional copula is rather given in terms of the inverse generator *ψ* −1 *y* (*s*) = *ψ* −1 *sH*(**y**) − *ψ* −1 *H*(**y**) , see (Charpentier 2003; Spreeuw 2006; Sungur 2002)) Applying (4) to the SAHARA family gives

$$\psi\_{y\cdot\theta,\varepsilon}(s) = \left(\frac{s+t-\varepsilon+\sqrt{\delta^2+\left(s+t-\varepsilon\right)^2}}{t-\varepsilon+\sqrt{\delta^2+\left(t-\varepsilon\right)^2}}\right)^{-(1+1/\theta)} \left(\frac{s+t-\varepsilon+\left(1+1/\theta\right)\sqrt{\delta^2+\left(s+t-\varepsilon\right)^2}}{t-\varepsilon+\left(1+1/\theta\right)\sqrt{\delta^2+\left(t-\varepsilon\right)^2}}\right),\tag{5}$$

so the conditional copula is again SAHARA with parameters *θ* and *e* − *t*. It follows that dependence increases over time, and the copula converges to Clayton with parameter *θ*. Again, this is unlike most other copula families where the limiting dependence is either none (independence) or perfect positive.

Some scatterplots follow in Figures 1–3, for Kendall's tau fixed at 0.25 and varying values for *e* and *θ*. As *e* went up, we encountered on the one hand increasingly positive dependence in the bottom left part, and increasingly negative dependence in the top right half. Such families could be considered when data feature strongly positive dependence for small values, and weakly positive, no, or even negative dependence for large values.

The SAHARA copula is clearly flexible and versatile. A drawback is that the inverse of the generator is not available in closed form, like some families introduced in McNeil and Nešlehová (2010), Hua and Joe (2011) and Hua (2015), rendering computations more complicated.

**Figure 1.** Scatterplot of the SAHARA copula for *θ* = 1 and *e* = 0; *τ* = 0.25.

**Figure 2.** Scatterplot of the SAHARA copula for *θ* = 4.5464 and *e* = 2; *τ* = 0.25.

**Figure 3.** Scatterplot of the SAHARA copula for *θ* = 69.11 and *e* = 10; *τ* = 0.25.

#### **4. Application**

For the numerical application in this section, we use the example about modelling dependence of coupled lives in Luciano et al. (2008) and Spreeuw (2014). The two publications used different data, although they both concern samples specified as generations from the same large dataset of annuitants from a Canadian insurer. In this section, copula families are fitted into the data from Luciano et al. (2008) rather than those from Spreeuw (2010). We follow the same procedure of modelling and calibration as that in Luciano et al. (2008) and Spreeuw (2014). Some elaboration on deriving the empirical generator is in order to render this paper self-contained.

The joint survival function of two remaining lifetimes *T m x* (male, age *x* at the start of the observation) and *T f <sup>y</sup>* (female, age *y* at the start of the observation) is given in terms of a survival copula *Cxy* as

$$\mathcal{S}\_{xy}(s,t) = \mathcal{C}\_{xy}(\mathcal{S}\_x^m(s), \mathcal{S}\_y^f(t)).$$

In this setup, the lives are coupled at the time when they are observed (rather than at birth, as in, e.g., Frees et al. (1996)), just like in Carriere (2000). Using a modified version of the procedure by Wang and Wells (2000), the performance of a candidate Archimedean copula is judged on the basis of distance between the empirical Kendall function, denoted by *K*b *n*(*xy*) , and the theoretical Kendall function, denoted by *K<sup>ψ</sup>* −1 **A**b (*xy*) (*v*), where *ψ* −1 **A**b is the inverse generator of the copula concerned, with **A**b being the parameter vector estimate minimizing the distance between *K*b *n*(*xy*) and *K<sup>ψ</sup>* −1 **A**b (*xy*) (*v*). For single parameter copulas, **A** = *θ*, while for families with two parameters, **A** = {*θ*, *e*}. The distance or error is defined under the *L* <sup>2</sup> norm (so, in the usual quadratic sense). Therefore,

$$\operatorname{error}\left(\psi\_{\widehat{\mathbf{A}}}^{-1}(\mathbf{x}y)\right) = \int\_{\mathfrak{f}}^{1} \left(K\_{\psi\_{\widehat{\mathbf{A}}}^{-1}(\mathbf{x}y)}(v) - \widehat{K}\_{n(\mathbf{x}y)}(v)\right)^{2} dv\_{n}$$

with

$$
\widehat{\mathbf{A}} = \arg\min\_{\mathbf{A}} \int\_{\xi}^{1} \left( \boldsymbol{K}\_{\boldsymbol{\psi}\_{\mathbf{A}}^{-1}(\mathbf{x}y)}(\boldsymbol{v}) - \widehat{\boldsymbol{K}}\_{n(\mathbf{x}y)}(\boldsymbol{v}) \right)^{2} d\boldsymbol{v}.
$$

Given that the data were right censored, the lower bound *ξ* was greater than zero. In this example, it is taken to be the smallest value for which *K*b *n*(*xy*) is positive:

$$\xi = \min \{ \nu : \widehat{\mathbb{K}}\_{n(xy)}(\upsilon) > 0 \}.$$

The empirical Kendall function, denoted by *K*b *n*(*xy*) , was derived from Dabrowska's nonparametric estimator of the joint survival function (see Dabrowska (1988)). Given that the data were right censored, with many observations being doubly censored, *K*b *n*(*xy*) is zero between 0 and a certain value *ξ*<sup>1</sup> > 0, at which point it jumps. In this case, *ξ*<sup>1</sup> = 0.23. The pseudo-maximum likelihood (PML) procedure uses as input rescaled Kaplan–Meier estimates of the marginal survival functions in order to accommodate censoring.

Luciano et al. (2008) fit the data to Clayton, Gumbel-Hougaard, Frank, entry 20 of Table 4.2 in Nelsen (2006) ("4.2.20 Nelsen") and the so-called Special copula. For convenience, the last two are listed below with their generators:

1. 4.2.20 Nelsen: *ψ<sup>θ</sup>* (*t*) = {log[*e* + *t*]} − 1 *<sup>θ</sup>* , *θ* > 0.

$$\text{2.} \quad \text{Special:} \quad \psi\_{\theta}(t) = \left(\frac{-t + \sqrt{t^2 + 4}}{2}\right)^{\frac{1}{\theta}}, \quad \theta > 0.$$

Luciano et al. (2008) concluded that 4.2.20 Nelsen fit the data best. In this section, we compare its performance with that of SAHARA and the best contender of common two-parameter families from Joe (1997, 2015), i.e., BB2. Its generator is:

$$\psi\_{\theta}(t) = \left\{ 1 + \frac{\log[1+t]}{\varepsilon} \right\}^{-\frac{1}{\theta}}, \quad \theta, \varepsilon > 0.$$

The 4.2.20 Nelsen is a special case of BB2 arising for *e* = 1.

Results are given in Table 1. The positive estimate for *e* indicates the absence of SI.

**Table 1.** Results for several copula families.


As in Luciano et al. (2008), we performed a graphical comparison between the theoretical and the empirical K functions through transformation *λ*(*w*) = *w* − *K*(*w*) for *ξ*<sup>1</sup> ≤ *w* ≤ 1. The result can be found in Figure 4. SAHARA achieved a significant improvement to the fit compared to the other families.

**Figure 4.** Graphical comparison between theoretical *λ*(*w*) = *w* − *K*(*w*) for SAHARA (yellow), BB2 (green) and 4.2.20 Nelsen (orange) and empirical one (blue).

Now consider the notion of SI in more detail. For two random variables *X*<sup>1</sup> and *X*2, *X*<sup>2</sup> being SI in *X*<sup>1</sup> is equivalent to Pr[*X*<sup>2</sup> > *x*2|*X*<sup>1</sup> = *x*<sup>1</sup> ] being nondecreasing in *x*<sup>1</sup> for all *x*2. Related to this is the notion of long-term dependence as introduced in Hougaard (2000). If we define *µ m t T f <sup>y</sup>* = *t<sup>y</sup>* as the conditional force of mortality of life (*x*) at duration *t* given *T f <sup>y</sup>* = *t<sup>y</sup>* < *t* ((*y*) dies at duration *ty*), then the dependence between *T m <sup>x</sup>* and *T f y* is of the long-term type if *µ m t T f <sup>y</sup>* = *t<sup>y</sup>* is constant or decreasing as a function of *ty*, while dependence is short-term if *µ m t T f <sup>y</sup>* = *t<sup>y</sup>* is increasing as a function of *ty*. To understand this, it is important to that, as indicated before, for Archimedean copulas, stochastic increasing (SI) is equivalent to −*ψ* 0 being logconvex. Spreeuw (2006) showed that this property of the generator also applies to long-term dependence, implying that SI and long-term dependence are actually equivalent. On the other hand, however, for the SAHARA family, we have:

$$\frac{\partial \left\{ \ln \left[ -\psi\_{\theta,\epsilon}'(s) \right] \right\}}{\partial s} = -\frac{1 + 1/\theta}{\sqrt{\left(s - \epsilon\right)^2 + 1}} \epsilon$$

which is monotone increasing in *s* across the board for negative *e*. For positive values of *e*, however, the expression decreases in *s* for 0 ≤ *s* < *e*, so the SI property does not hold. The positive parameter estimate for *e* suggests that short-term rather than long-term dependence may prevail between the coupled lives. To investigate this further, we analysed the data in the same vein as in Spreeuw and Owadally (2013), who devise an augmented Markov model to allow for short-term dependence for the entire dataset. Results are reported in Table 2.

In this table, *e* denotes the time in which an integer number of years that have elapsed since the death of the partner. So, e.g., *e* = 0 concerns the lives that were bereaved less than a year ago. For each possible value of *e* (noting that each life was observed for 5 years or less), we calculate number of deaths reported, the risk exposure, and the overall mortality rate being the ratio of the values in the second and third column. So, for instance, the risk exposure of lives who lost their partner less than a year ago is 604.87, and there were 69 lives that died within one year after their partner. Now, long-term dependence implies that the mortality rate in the last column should be increasing as a function of *e*, but the results in Table 2 show that this is not the case, and that short-term dependence may be present even though the aggregate mortality rate for *e* = 4 is higher than for *e* equal to 1, 2 or 3.

Ideally, one such table should be shown for each gender. However, due to the small number of observed deaths in the dataset, in particular for higher values of *e* caused by heavy censoring, males and females were combined. Results in Table 2 should thus be interpreted as an indication of possible short-term dependence, rather than firm evidence.


**Table 2.** Mortality for all couples, with *e* denoting the number of years since partner's death.

#### **5. Conclusions**

In this paper, we introduced a new Archimedean copula family derived from the SAHARA utility function. With SAHARA utility first increasing to a maximum and subsequently decreasing, the corresponding copula family allows for stochastically increasing (SI) and non-SI at the same time, depending on the sign of one of the parameters. As the numerical application shows, this family could fit the mortality data of coupled lives well. The parameter estimates suggest the possible existence of short-term dependence, i.e., the mortality of bereaved lives increases on bereavement but diminishes later.

**Funding:** This research received no external funding.

**Data Availability Statement:** Anyone who would like to use the data should consult Professor Edward W. (Jed) Frees of University of Wisconsin-Madison, USA, e-mail: jfrees@bus.wisc.edu, or Professor Emiliano Valdez of University of Connecticut, USA, e-mail: emiliano.valdez@uconn.edu.

**Acknowledgments:** We wish to acknowledge the Society of Actuaries, through the courtesy of Professors Edward (Jed) Frees and Emiliano Valdez, for providing the data used in this paper.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Risks* Editorial Office E-mail: risks@mdpi.com www.mdpi.com/journal/risks

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-8390-7