Global Optimization for Automatic Model Points Selection in Life Insurance Portfolios

: Starting from an original portfolio of life insurance policies, in this article we propose a methodology to select model points portfolios that reproduce the original one, preserving its market risk under a certain measure. In order to achieve this goal, we ﬁrst deﬁne an appropriate risk functional that measures the market risk associated to the interest rates evolution. Although other alternative interest rate models could be considered, we have chosen the LIBOR (London Interbank Offered Rate) market model. Once we have selected the proper risk functional, the problem of ﬁnding the model points of the replicating portfolio is formulated as a problem of minimizing the distance between the original and the target model points portfolios, under the measure given by the proposed risk functional. In this way, a high-dimensional global optimization problem arises and a suitable hybrid global optimization algorithm is proposed for the efﬁcient solution of this problem. Some examples illustrate the performance of a parallel multi-CPU implementation for the evaluation of the risk functional, as well as the efﬁciency of the hybrid Basin Hopping optimization algorithm to obtain the model points portfolio.


Introduction
This work deals with the risk management in life insurance companies, which is a key aspect in the insurance industry and in Solventia II regulation. Asset Liability Management (ALM) is the process of managing the assets and liabilities portfolios as well as cash flows to reduce the firm's risk of loss. It plays a key role in the profitability of insurance companies, see [1,2] and the references therein. From the numerical point of view it mainly consists of computing the joint projections of the future cashflows of assets and liabilities portfolios, which can be done by using Monte Carlo algorithms.
These kind of portfolios comprise a huge number of policies with several characteristics: maturity, possibility of early cancelation, age of the policy holders, gender, etc. Computing this stochastic projection of the balance sheet related to these portfolios, involves using large number of scenarios with the corresponding interest rate models, mortality model, etc. This projection can be computed using the whole number of policies, what is referred as stand alone projection, where each policy is projected in each scenario. In practice, the computation of the projections of the liabilities in original portfolios (containing hundreds of thousands of policies) often becomes a highly demanding computational time objective, making the problem barely tractable even if High Performance Computing (HPC) hardware is employed.
For that reason, insurance companies perform these projections by grouping similar policies in some representative contracts, known as model points. The choice of this grouping of policies must preserve the representation of the risk distribution of the original portfolio. We refer to the document [3] for further details.
Although model points have been used by life insurance companies for a long time, the more theoretical mathematical foundations are not so developed to our knowledge. Particularly, a rigorous framework to measure the risk associated to the replacement of the original portfolio by the model points portfolio is an important task that has not been much addressed in the literature. For example, in the interesting work of [1,2] about ALM for an insurance company, model points portfolios are just defined without being built from original policies portfolios. More recently, in [4], a conservative model points portfolio is performed by controlling the Tail-Var risk measure.
Note that the new model points portfolio provides a highly reduced representation of the original portfolio that should preserve the same risk properties. For this purpose, an appropriate function to measure the associated risk in this replacement has to be chosen.
In the present work, the model points portfolio selection problem is framed in the domain of portfolio representation, where the risk functional gauges the averaged differences between the values of original and model points portfolio, the last one taken in a specific set of portfolios [5]. This is one of the main innovative parts of the present work. The risk measure is defined attending to the variation of some of the stochastic underlying risk components for a certain time horizon. The models for the interest rate evolution and the population mortality become two of the key factors in the risk management of life insurance portfolios [6]. In the present paper, we just consider the uncertain future fluctuations of interest rates. Although other alternative models could be considered, such as the popular short rate models, we have chosen a classical version of the LIBOR (London Interbank Offered Rate) market model for forward LIBOR rates. Additional risk factors could be incorporated in future works, as for example mortality risk, as well as the dependence among them when more than one underlying risk factor is considered.
However, even after choosing the more appropriate risk measure functional, finding the structure of the model points portfolio regarding this risk measure can turn out to be a extremely hard computational task. As we show in this article, this problem can be written a global optimization one, consisting in the minimization of the distance between the real and the model points portfolios, once the appropriate risk measure has been chosen.
The main goal of this article is to develop an original technique for the automatic generation of the model points portfolio. This objective mainly involves two main tasks. First, we have to choose/fix and build a correct risk measure attending to market models. In our case we will build this measure by using the LIBOR market model for forward interest rates (see [7], for details). Secondly, once we have formulated the problem in terms of a distance minimization based on this measure, we need to numerically compute this objective function in a highly efficient way, as next we have to build a numerical optimization method for the minimization of this objective function. As we will see, the resulting optimization problem is a global optimization one, so that efficient global optimization algorithms are required. More precisely, the algorithms need to be fast and robust, so that they won't get stuck in local minima.
Moreover, we note that for large real life policies portfolios, the here proposed methodology involves a huge computational cost. In this case, specific HPC techniques should be used. In [8] a full HPC implementation is carried out by combining different parallelization techniques for multi-core and many-core architectures.
The plan of the paper is as follows. In Section 2, we briefly summarize the main issues related to the LIBOR market model that is one of the main point in the definition of the risk functional. In Section 3, we formulate the problem of the model points portfolio selection. In Section 4, we define the risk measure here proposed. In Section 5, we describe the discretization of the cost function. Section 6 contains some numerical tests to validate the convergence of the optimization algorithm to the analytical solutions and finally we show an application to a real portfolio without analytical solution.

Model Setup LIBOR Market Model
As we mentioned in the introduction, we just consider the interest rate as the only risk factor. In this section, following [7] we introduce the dynamics of the discounted bond price, as it follows from the classical LIBOR Market Model governing the forward rates evolution in time. It is important to emphasize that any of the existing stochastic models for the evolution of interest rates could be considered in our proposed methodology and its corresponding Monte Carlo simulation can be included in the developed software toolbox as a building block of the global algorithm, that can take advantage of the proposed parallelization as well. Thus, classical short rate models as Vasicek, CIR or Hull-White, as well as models incorporating possible negative rates as shifted-SABR, free boundary SABR or displaced diffusion can be considered. The volatility smile could be incorporated through the more complex SABR-LIBOR models, for example. The choice of LIBOR Market Model provides an additional computational complexity and cost in the simulation of the risk factor which enters in the evaluation of the risk functional. Although we are aware of the evolution of this model and the difficulties to cope with the practical presence of a multicurve setting (see [9][10][11], for example), we have chosen its classical version as described in [7].
Thus, in this section we first summarize the main concepts of the classical LIBOR Market Model [7] to be used in the sequel. Next, we introduce the discounted bonds price, for which we prove the dynamics satisfied by the process in the final lemma of the section.
Libor rates. For a given fixed integer N > 0 we define the set T = {T 0 , T 1 , ..., T N } that represents the fixed tenor structure, with T 0 = 1 and T 0 < T 1 < ... < T N , being T n the maturity. We denote τ n = T n − T n−1 , for n = 1, ..., N, the corresponding accruals. Moreover, we set the one year time interval I = (0, 1).
Here and in the sequel, we denote by B n (t) the discounted price of a zero-coupon bond at time t ∈ I under the risk neutral measure. This bond matures at the tenor date T n , for n = 0, ..., N. Moreover, let F n (t) be the value at time t ∈ I of the (simply-compounded) LIBOR forward rate that operates during the accrual period (T n−1 , T n ], for n = 1, ..., N. We recall that F n (t) is defined in such a way that the following condition is met: Hence, by recursively applying the previous condition, we have for any t ∈ I.
Since t < T n , then the price B n (t) is given by the previous expression, for any t ∈ I and n = 0, ..., N.
Stochastic setting. Hereafter we consider B 0 (t) the price of the bond with expiring tenor date T 0 = 1, for t ∈ I, as the process acting as reference numeraire. Furthermore, let Q be the forward measure corresponding to T 0 , that is, the martingale measure associated to the numeraire B 0 .
Moreover, we shall fix a N-dimensional Wiener process W(t) = (W 1 (t), ..., W N (t)), for t ∈ I, defined on the suitable complete probability space (Ω, F , Q), and we introduce the associated correlation matrix = ( nk ) nk , so that dW n (t)dW k (t) = nk dt.
More precisely, we consider constant correlation coefficients between the LIBOR forward rates, given by the parameterization: Moreover, for any h = (h 1 , ..., h N ) ∈ R, we define the following norm: Any stochastic process shall be defined on the probability space (Ω, F , Q). Further, since Q is understood as the reference market measure, we shall refer to any process defined on such a probability space as a risk-free process.
LIBOR Market Model. For each n = 1, ..., N, let σ n be a given deterministic function defined on the interval I that represents the volatility of F n . Given any fixed n = 1, ..., N, the LIBOR Market Model assumes that the risk-free dynamics associated to the process F n (t), for t ∈ I, satisfies dF n (t) = µ n (t)dt + σ n (t)F n (t)dW n (t), jointly with a given initial value of the forward rate, F n (0), where the drift component µ n (t) at any time t ∈ I is completely determined by following identity: We choose the usual parameterization for the volatilities σ n (t): For any fixed t ∈ I, we set µ(t) = (µ 1 (t), ..., µ N (t)) and define Σ(t) as the matrix: Σ nk (t) = σ n (t)F n (t)δ nk , for any n, k = 1, ..., N, being δ nk the Kronecker delta. Besides, we use Σ n (t) to denote the n-th row of Σ(t), for n = 1, ..., N. Then, when setting F(t) = (F 1 (t), ..., F N (t)), we take (1) as a N-dimensional dynamics using the following compact notation: jointly with the initial condition F(0) = (F 1 (0), ..., F N (0)). Moreover, the discounted price process linked to the bond expiring at tenor date T n , is denoted byB for any t ∈ I, for any n = 1, ..., N.
The dynamics under the risk neutral measure of the discounted price for any bond that expires at the tenor date T , is given in the following result. Lemma 1. For any n = 1, ..., N, the dynamics of the discounted bond price processB n (t) satisfies the following stochastic differential equation: with ε n (t) given by Proof. Fix n = 1, ..., N and notice that according to the identity (1) the discounted pricẽ p n (t) at time t ∈ I, that is given by (4), satisfies the condition For any diffusion process X(t), for t ∈ I, driven by the Wiener process W, let DCX(t) denote its vector diffusion coefficient. In particular, the representation (3) implies that DCF n (t) = Σ n (t), (6) for any t ∈ I. Thus, we obtain Finally, taking into account that for any t ∈ I, the proof is completed.

Insurance Policies Portfolio
Life insurance contracts usually can pay the insurer either a stream of cash flows during the lifetime of the policyholder or a unique lump sum benefit that is paid upon her/his death under certain conditions. When agreed in the clauses of the contract, other situations as the appearance of a terminal illness might also lead to early payments. The policy is alive while the policyholder pays a specific premium, that may be either in the form of periodic payments or consist of an initial lump sum. We assume to deal with a life insurance policy that is not affected by credit risk, which means that the insurance company always guarantees the entire benefit that is settled in the contract.
In this section, we describe the model points selection problem for the case of a portfolio of term insurance policies. A term insurance policy pays a lump sum benefit fixed in the contract, in the case that the policy owner death happens before an specific term fixed in the contract. In order to simplify the presentation, the benefit associated to each policy is assumed to be represented by a unit amount of a given currency, and that the policies are not affected by credit risk. Furthermore we do not consider the premiums, nor any further expenses that are responsibility of the client.
Thus, in this section, after briefly describing the setup for the insurance portfolio we consider, we mainly introduce two new definitions that allow to rigorously pose the problem in the mathematical framework.
Term insurance portfolios. We assume that each term insurance policy inside a portfolio is identified by the policy owner age at time t = 0 and the term date of the contract. In this respect, let us fix I, J ∈ N and let X = {x 1 , ..., x I } and Y = {y 1 , ..., y J } be two finite sets of real values such that x i ≥ 0 and y j ≥ 1, for any i = 1, ..., I and j = 1, ..., J. Here and in the sequel, any couple (x i , y j ) uniquely defines the family of policies related to the class of individuals that are aged x i ∈ X at time t = 0 and to the term of the contract y j ∈ Y. Thus, the original insurance portfolio is a set of I × J policies.
Next, for any i = 1, ..., I, we denote by µ(s, x i + s) the mortality force at time s ≥ 0 associated to the individuals with age x i ∈ X . Moreover, we suppose that the function µ(s, x i + s) is deterministic and observable, taking values for any x i ∈ X and s ≥ 0. Furthermore, the following convenient condition is assumed: for any i = 1, ..., I and any s ∈ I. (8) Remark 1. Note that hypothesis (8) guarantees that any portfolio of term insurance policies remains unchanged within the time interval I due to the policy owners death. Besides, hypothesis (8) is acceptable since the performance of the overall portfolio is barely affected by the events taking place during the first year.
According to this setup, we define the survival index as which provides the proportion of persons aged x i that survive up to age x i + T n . We consider the Gompertz-type law modelling the force of mortality [12]: where a and b are given observable functions defined for s ≥ 1. Assuming that a(s) = 0 for s ∈ I we guarantee that the previous convenient assumption (8) holds. In what follows, we use S T to denote the derivative of S with respect to its second variable: Now, we introduce the definition and notations for term insurance policies.

Definition 1.
For any x i ∈ X and y j ∈ Y, the discounted risk-free value at time t ∈ I of a term insurance policy owned by an individual with age x i ∈ X and with term y j ∈ Y is given by Any linear combination of the processes (10) gives the risk-free discounted value of a term insurance portfolio. This is stated by the following definition. Definition 2. We call (term insurance) policy portfolio any matrix v = (v ij ) ij , for i = 1, ..., I and j = 1, ..., J. Moreover, for any policy portfolio v, its risk-free discounted value v(t) at time t ∈ I is defined as Given any policy portfolio v, we understand any of its components v ij as the amount of policies owned by the class of individuals labeled by the age x i ∈ X and the term y j ∈ Y. In what follows, for any couple of policy portfolios v 1 and v 2 we shall write , for any t ∈ I.

Model Points Risk Estimation
In this section, we present one of the main innovative issues of the present article, which mainly consists in the application of the theory presented in [5] to the previously introduced LIBOR rates setting and the portfolio of term insurance policies. More precisely, we take into account LIBOR forward rates dynamics as the stochastic underlying stochastic factor. For this purpose, we shall always assume a policy portfolio v relative to X and Y to be a priori given. Moreover, we consider a fixed set of term insurance policy portfolios W that contains policies with L given ages and M given maturities, so that L is smaller than I and M is smaller than J. Then, each w = (w lm ) ∈ W is referred as one model points portfolio in W. In this setting we introduce the following definition. Definition 3. Let R(·|v) : W → R be the model points risk functional induced by v over W, which is defined as the follows for any w ∈ W, We understand the model points risk functional R(w|v) as a measure of the error when representing (the original policies portfolio) v by (the model points portfolio) w ∈ W. This error is computed as the average of the difference between both portfolios along time, in terms of the variation of the stochastic underlying factors to be considered, which is the interest rate risk in our case.

Definition 4.
Let R(·|v) be the model points risk functional induced by v over W. A model points portfolio w * ∈ W is defined to be optimal with respect to v if the following inequality holds true for any w ∈ W. (12) Notice that, in the special case when v ∈ W one has that R(v|v) = 0. In this respect, we regard any model points portfolio w * ∈ W satisfying the inequality (12) as an optimal representation of the policy portfolio v.
The following result provides an alternative representation of the model points risk functional induced by v over W.
Proposition 1. The model points risk functional R(·|v) is given by: Proof. Let E = R N and letB(t) = (B 1 (t), ...,B N (t)) for any t ∈ I. Moreover, let U = (X × Y ) R be the class of real matrices r = {r(x i , y j ) : x i ∈ X , and y j ∈ Y }.
Consider the functional Z ∈ L(E, U) which is given for any q ∈ E by and note that the process z(t) = {z(t) : t ∈ I} defined by (10) is obtained for the choice , for any t ∈ I.
Note that Z(q) is linear, for any q ∈ E. Then, the Frechét derivative denoted by ∇Z(q) ∈ L(E, U) satisfies the following identity: for x i ∈ X and y j ∈ Y.
for any q ∈ E. On the other hand, for any t ∈ I, when letting ε(t) be the matrix with the n-th row given byB n (t)ε n (t), where ε n (t) is defined in (5), we have that for any x i ∈ X and y j ∈ Y, the following identity holds Consider now the function ζ : I × E → U, which is defined by for any t ∈ I and q ∈ E.
As condition (8) is imposed, the survival index (9) is independent on t ∈ I, so that the representation (15) is allowed.
Then, since (4) may be rewritten as the following E-valued dynamics, then Proposition 5 in [5] can be applied to obtain for any w = (w lm ) ∈ W. (17) which, jointly with the identity (14), gives the representation (13).
Therefore, the construction of a model points policies portfolio w * for an original portfolio v is given by: As a consequence, obtaining the model points portfolio is a very hard problem, as it is given by the solution of a global optimization problem in high dimension. Note that the number of policies in the model points portfolio, I × J, defines the dimension of the problem. Therefore, efficient numerical methods are required to minimize the model points risk functional R(· | v) given by (13). Note that this risk functional is referred as the cost function in the literature related to optimization methods.

Numerical Methods
In this section, we describe the two main involved tasks in the numerical methods, which represent innovative and original contributions of the present article. On the one hand we introduce the appropriate discretization of the model points risk functional (13) and its parallel implementation. On the other hand, we choose the global optimization numerical methods for solving problem (18), which require a very efficient evaluation of the model points risk functional. For this purpose, a parallel multi-CPU implementation for the evaluation of the risk functional has been developed from scratch. This implementation runs inside the hybrid Basin Hopping optimization algorithm to select the model points portfolio. To our knowledge, this is the first time these HPC computational tools are applied to this kind of problems to speed up the computational time.
The expectation that appears in the risk functional (13), is obtained by means of Monte Carlo techniques. Note that each simulation involves the computation of forward LIBOR rates evolution that follow the previously described LIBOR Market model. Following the ideas in [7], taking logarithms and using Ito lemma in (1), the forward rate dynamics is given by the equation: with a deterministic diffusion coefficient. Now, for the time discretization of (19), we apply the Euler-Maruyama numerical method, thus obtaining for n = 1, . . . , N, whereF n (t) is the approximation to F n (t) andŴ n (t + ∆t) −Ŵ n (t) ≡ √ ∆tN (0, 1) corresponds to the discretizaton of dW n (t). We consider a uniform mesh with nodes t q = q∆t, for q = 0, . . . , N t , with ∆t the constant time step.
The discretization of expression (13) for the functional,R(w|v), is given by: with C representing the correlation matrix, N t being the number of time steps and N p denoting the number of Monte Carlo simulations. Moreover, the expression of R pq is as follows: with the vector ν pq (T n ) given by In the previous expression the index p corresponds to a particular simulation of the forward LIBOR rates, index q is associated to time t q and n is linked to the maturity T n . Furthermore, let us consider the notation v ij being the nominal of policies corresponding to age x i and to maturity y j in the initial portfolio, while w lm represents the same quantity in the portfolio of model points.
Therefore, the computational cost associated to the cost function (given by the risk functional (13)) evaluation turns to be actually high: note that we use a Monte Carlo numerical scheme with each path involving the simulation of forward LIBOR rates and additionally we have a computationally intensive loop in the policies corresponding to the original portfolio.
Thus the cost function needs to be calculated in a very efficient way, because any global optimization algorithm essentially requires the repeated evaluation of the cost function. For this purpose, we have carried out the parallel implementation of the risk functional (13), by using OpenMP in a multi CPU architecture. A parallel random number generation algorithm has been used for the parallelization of the Monte Carlo algorithm.
More precisely, in our case we have used the well known Tina's Random Number Generator library (TRNG). The C++ code for the LIBOR interest rates generator, together with the Monte Carlo algorithm for the risk functional, has been written from scratch.
For solving the global optimization problem, we use stochastic algorithms. One example of these kind of algorithms is Simulated Annealing (SA): this algorithm consists in exploring the search space by generating a sequence of neighbor points, where the cost function is evaluated. The points are accepted or discarded according to a physical based acceptance criterion, like Lattice-Boltzmann's law, see [13] and references therein for details. Global metaheuristic pure stochastic algorithms, like SA, have slow convergence and thus a higher computational cost. One alternative is the use of deterministic local optimization algorithms, but although they are faster, they could be trapped into a local minimum. In this work, we consider also hybrid algorithms, to benefit from the speed of local deterministic algorithms, while retaining the global convergence of stochastic ones. More precisely, we use the Basin Hopping (BH) algorithm, see [14,15], jointly with the references they include. Basically, inside the BH algorithm, a Simulated Annealing is applied for randomly sampling the search space, while local gradient optimizers are launched from the SA randomly generated points. In this way the convergence of SA is accelerated. We use L-BFGS-B as the local optimizer in BH. The whole optimization routines have been implemented from scratch in C++, following the previous works [13,15].

Numerical Experiments
In this section, first we present two tests with known solution to validate the proposed model points selection technique, and finally we show an application to a real world scenario. For each experiment we show the comparison, in number of function evaluations, between SA and BH.
For all the experiments we use the same settings for the LIBOR interest rate model and the Monte Carlo numerical method. Concerning the LIBOR model: we take β = 0.01 for the correlation coefficients between the forward rates; and we take 100 tenors with maturities between 1 and 100, and initial rates F 1 (0) = 0.01, F 2 (0) = 0.02, F 3 (0) = 0.03, F 4 (0) = 0.04 and F n (0) = 0.05, for n ≥ 5. We consider the values a = 0.07, b = 0.2, c = 0.6 and d = 0.075 for the parameterization of the volatilities. For the mortality force, we consider constant parameters, µ(x) = a exp(bx), where we will take a = 0.0003 and b = 0.06. From the numerical point of view, all the tests in this section have been performed by using 1000 paths for the Monte Carlo simulation for the evaluation of the cost function.
Regarding to the hardware specifications, all tests have been carried out in a server with 16 CPU cores (two eight-core Intel Xeon E5-2620 v4).

One Model Point Portfolio
The goal of this example is to represent an original portfolio with a model points portfolio with one single model point.
More precisely, the original portfolio contains 10 policies, see Table 1. The model points portfolio consists in just one single model point with age 42 and maturity 28 years. The objective is to determine the total nominal in this model point that minimizes the corresponding risk functional (13).  Therefore, in this case the cost function (risk functional) is given by a one dimensional function. If we evaluate the cost function for different number of contracts, we obtain the expected parabolic convex function, so that the global minimum can be graphically checked (see Figure 1). Therefore, in this case the cost function (risk functional) is given by a one dimensional function. If we evaluate the cost function for different number of contracts, we obtain the expected parabolic convex function, so that the global minimum can be graphically checked (see Figure 1).

Repeated Policies
In this test we use a problem with known solution to assess the proposed methodology. More precisely, we check and compare the convergence of SA and BH algorithms to the exact solution. Furthermore, we show the speedup of the parallel version of the cost function.
To build the original portfolio, we repeat each policy in the portfolio of Table 1 1000 times, so that we obtain a portfolio containing 10, 000 policies. Although they are repeated, each policy is treated as a different individual one. The goal is to represent the original portfolio with the 10 model points portfolio in Table 1. Therefore, in this problem we aim to find the nominals in the model points portfolio that better mimic the original portfolio, regarding the risk functional. In this case the analytical expression for the solution can be trivially found by multiplying the number of contracts in Table 1 by 1000 and the corresponding value of cost function is equal to zero.
As we can see in Figure 2 for the BH method with L-BFGS, the method reaches the very small value 1.75 × 10 −8 of the cost function. In Table 2 we show the obtained values of the nominals for the model points portfolio rounded to the eighth decimal digit. The resulting difference between 1.75 × 10 −8 and zero in the cost function is due to this rounding procedure. The computational time was 87.35 h (about 3.63 days) using L-BFGS, when using 16 cores (and 32 threads).

Repeated Policies
In this test we use a problem with known solution to assess the proposed methodology. More precisely, we check and compare the convergence of SA and BH algorithms to the exact solution. Furthermore, we show the speedup of the parallel version of the cost function.
To build the original portfolio, we repeat each policy in the portfolio of Table 1 1000 times, so that we obtain a portfolio containing 10,000 policies. Although they are repeated, each policy is treated as a different individual one. The goal is to represent the original portfolio with the 10 model points portfolio in Table 1. Therefore, in this problem we aim to find the nominals in the model points portfolio that better mimic the original portfolio, regarding the risk functional. In this case the analytical expression for the solution can be trivially found by multiplying the number of contracts in Table 1 by 1000 and the corresponding value of cost function is equal to zero.
As we can see in Figure 2 for the BH method with L-BFGS, the method reaches the very small value 1.75 × 10 −8 of the cost function. In Table 2 we show the obtained values of the nominals for the model points portfolio rounded to the eighth decimal digit. The resulting difference between 1.75 × 10 −8 and zero in the cost function is due to this rounding procedure. The computational time was 87.35 h (about 3.63 days) using L-BFGS, when using 16 cores (and 32 threads).  Next, in Figure 3 we show the convergence of the SA algorithm for solving the optimization problem. In this case, the computational time was 839.35 h (about 35 days). The obtained solution is the same as with BH algorithm.
In this test, we have validated that the proposed methodology can group the repeated policies in the corresponding bucket, which is an expected property of the risk functional. In view of computational times we note that BH results to be more efficient than SA. We also have shown that the problem has high computational cost.
We can also use this test to show the performance of the multi-CPU implementation of the model points risk functional calculation. For this purpose, we consider from 100 to 10, 000 policies in the original portfolio, by repeating the policies of the portfolio in Table  1. Moreover, we consider 1000 paths for the Monte Carlo simulation. Figure 4 and Table  3 illustrate the achieved speed up with the parallel implementation of the risk functional with different numbers of policies in the original portfolio and different number of cores being used. As it is illustrated by both, it occurs that the speed-up gets closer to the number of used cores when the number of policies increases. Clearly, the larger the portfolios the more advantageous and efficient the parallelization turns out to be.   Next, in Figure 3 we show the convergence of the SA algorithm for solving the optimization problem. In this case, the computational time was 839.35 h (about 35 days). The obtained solution is the same as with BH algorithm.
In this test, we have validated that the proposed methodology can group the repeated policies in the corresponding bucket, which is an expected property of the risk functional. In view of computational times we note that BH results to be more efficient than SA. We also have shown that the problem has high computational cost.
We can also use this test to show the performance of the multi-CPU implementation of the model points risk functional calculation. For this purpose, we consider from 100 to 10, 000 policies in the original portfolio, by repeating the policies of the portfolio in Table 1. Moreover, we consider 1000 paths for the Monte Carlo simulation. Figure 4 and Table 3 illustrate the achieved speed up with the parallel implementation of the risk functional with different numbers of policies in the original portfolio and different number of cores being used. As it is illustrated by both, it occurs that the speed-up gets closer to the number of used cores when the number of policies increases. Clearly, the larger the portfolios the more advantageous and efficient the parallelization turns out to be.

Realistic Portfolio
In this test, we consider a realistic synthetic portfolio that contains 10,000 different policies. In Figure 5 we show the ages, nominals and maturities of the first 20 policies of this portfolio, as showing all policies would require a very large table. Moreover, in Figure 6 we show the cashflows of the original portfolio due to the payments to the policyholders.   The target model points portfolio consists of 45 policies, that comprises a table with 9 ages and 5 maturities. Thus, the ages range from 30 to 70 years with a constant step of 5 years, while maturities range from 5 to 25 years with a constant step of 5 years.
Note that this test is much more difficult than the one in Example 2. The required total computing time by the BH hybrid optimization algorithm, with L-BFGS-B as local minimizer, was 136.43 h when using 16 processors for the parallel evaluation of the risk functional. Figure 7 illustrates convergence of the algorithm. After solving the optimization problem, the model points risk functional achieved the value 0.175927. Moreover, the computed nominals in the model points portfolio are first shown in Table 4 and Figure 8. Also, in Figure 9 we detail the same results in different graphics for each maturity. In the proposed hybrid algorithm, the global property becomes of great importance for this test, as several iterations of the global optimization algorithm were needed to reach the minimum, thus avoiding to get stuck to any of the possibly existing local minima as it could be the case if a pure local optimization method were applied. Also the convergence speed and accuracy of the L-BFGS-B local optimizer results to be a key point, as the SA pure global algorithm was not able to converge for this realistic example.

Conclusions
ALM is a key process to evaluate the profitability of insurance companies and to manage the assets and liabilities portfolios as well as cash flows with the aim to minimize the risk of loss of the company. In this ALM process, the computation of joint projections of future cashflows associated to both portfolios is mandatory.
On the liabilities side, portfolios contain a very large number of policies with different characteristics and the computation of the corresponding stochastic projection requires the use of large number of scenarios that involve the underlying stochastic factors (interest rates, mortality index, etc). Therefore, the projection of each policy in each scenario (stand alone projection) leads to highly demanding or prohibitive computational times, even when using HPC tools. In order to overcome this drawback, insurance companies group similar policies in some representative ones, known as model points. However, the model points portfolio must preserve the risk distribution with respect to the underlying factors of the original portfolio. Although model points have been used in practice, a rigorous framework to measure the risk associated to the replacement of the original portfolio by the model points portfolio has not been enough addressed in the literature. From the industrial side, this study has been jointly proposed in collaboration with the financial consulting company Analistas Financieros Internacionales (AFI) interested in the topic inside the EU H2020 MSCA-ITN-EID Wakeupcall project.
In view of previous summary of the state of the art which constitutes the first working hypothesis, we have mainly proposed a theoretically based methodology to address the selection of a model points portfolio for a given initial portfolio of life insurance policies. The methodology requires to overcome both mathematical and computational challenges.
A second hypothesis is that the model points portfolio selection requires definition of an appropriate market risk functional to evaluate how the it represents the distribution of risks associated to the initial one. Also the rigorous mathematical setting has to be posed. Assuming that the interest rates are the main underlying risk factor, an appropriate market risk functional has been proposed and the a rigorous mathematical framework coming from portfolio representation has been used. Although the LIBOR market model for the evolution of forward rates has been chosen, we point out that any other interest rate model can be plugged into our methodology, so that the choice of LIBOR market model is by no means restrictive.
Once the appropriate risk functional has been defined, obtaining the model points portfolio involves the solution of a high dimensional global optimization problem, so that the numerical optimization methods require a very large number of risk functional evaluations. The risk function has been discretized by using Monte Carlo technique. In order to speed up the high computational cost associated to each evaluation, it has been parallelized in a multi CPUs setting using OpenMP. Concerning the choice of the optimization methods, the hybrid Basin Hopping optimization algorithm, which combines the advantages of mixing stochastic global optimization algorithms with gradient local optimizers, has been proposed and compared with an alternative purely global optimization technique. The accuracy and performance of the algorithm have been analyzed with several numerical examples, both in cases with known solution and also in a real example without analytical solutions. The results show that the hybrid global optimization algorithm is the more suitable for selecting the model portfolio.
We have also shown that the developed methodology involves a huge computational cost for realistic large life policies portfolios, because the large number of policies and the highly demanding Monte Carlo method at each evaluation of the risk functional. This disadvantage could be overcome by using specific HPC techniques and architectures based on GPUs, see for example [8,16].
Although the main achievements are related to the proposed rigorous mathematical framework and methodology, as well as to the efficient computational implementation using HPC tools and parallelization of Monte Carlo technique for the risk functional evaluation, we also consider the possible interest of the results in the insurance sector. We note that the problem has been jointly proposed and assessed by professionals of the financial company AFI inside a the European Industrial Doctorate (EID) Wakeupcall project in the MSCA-ITN-EID call of H2020 Program. The obtained results provide a solid mathematical basis to the proposed methodology for building model points portfolios in life insurance policies and the efficient implementation allows to speed up the highly demanding computational cost, both objectives being of great advantage when addressing ALM in an insurance company.
Finally, we note that the optimal model points selection for life insurance portfolios has been addressed by considering that stochastic evolution of interest rates is the main source of risk. This means that no stochastic fluctuation affecting the time evolution of the mortality rate has been taken into account, as a deterministic survival trend over time for any cohorts is considered. However, evidences of the stochastic behavior in the trend of the life expectancy can be found in [17] and the references therein. In this respect, one may also refer either to mortality or longevity risk to contribute to losses and risks due to unexpected changes in the long-terms life trend displayed by the policy owners. As a direct consequence, a more sophisticated approach for the model points selection should be obtained by considering the stochastic time evolution of both the interest rate term structure and the mortality rate. In this context, the dependence between mortality risk and interest rate risk cover a central issue. A first attempt to model this dependence is carried out in [6] and future research could be done to incorporate this additional source of risk.