Several authors have investigated the insurance-linked derivatives valuation. The method usually employed is the development of pricing models based on the hypothesis of geometric Brownian motion, to systematize the instantaneous reporting claims evolution, and to incorporate the possibility of major catastrophes occurring through Poisson processes.

Cummins and Geman (1995) pricing of the first generation of catastrophe pricing of the first generation of catastrophe, Cat futures and Cat options, traded at the Chicago Board of Trade. They defined the loss index as the sum of the claims associated with each catastrophe and used a geometric Brownian motion to represents the randomness of the claims reporting process, and a Poisson process that incorporates the jumps in the claim process due to the occurrence of new catastrophes. Geman and Yor (1997), follow a similar approach but use the diffusion process with jumps to directly model the loss rate of Property Claim Services (PCS) options. Aase (1999) develop a valuation model of catastrophe futures when the loss index follows a stochastic process containing jumps of random claim sizes at random time points of accident occurrence. This model is a particular case of the model created by Embrechts and Meister (1997) which represents the behavior of the catastrophic loss index through a mixture of compound Poisson processes and a random loss frequency. Baryshnikov, Mayo and Taylor (2001), using continuous trading and risk neutrality, price the catastrophe bond using a double compound Poisson process to capture the different characteristics of catastrophe dynamics. Burnecki and Kukla (2003) apply the results of Baryshnikov, Mayo and Taylor (2001) to calculate non-arbitrage prices of a zero-coupon and coupon CAT bond. Muermann (2003) introduces the concept of actuarial consistency and derives a representation of the prices of non-arbitrage catastrophe derivatives (Cat futures and Cat options) written on an underlying loss index that is modeled as a compound Poisson process. Loubergé, Kellezi and Gilli (1999) pricing a loss index triggered cat bond applying the catastrophe option pricing model developed by Cummins and Geman (1995). Lee and Yu (2002) develop a contingent claim model to price a CAT bond through geometric Brownian motion. This model incorporates stochastic interest rates and considers moral hazard, basis risk and default risk. Biagini, Bregman and Meyer-Brandis (2008) value catastrophe options and describe the index using an inhomogeneous compound Poisson process for the loss period and use an inhomogeneous exponential Levy process to re-estimate the index during the development period and up to maturity.

Jaimungal and Wang (2006) analyze the pricing of catastrophic put options under stochastic interest rates with losses generated by a compound Poisson process. Asset prices are modeled through a jump-diffusion process which is correlated to the loss process. To evaluate these catastrophe options Wang (2016) employ a compound doubly stochastic Poisson process with lognormal intensity to describe accumulated losses and assume the volatility varies stochastically. Jarrow (2010) use a pricing methodology based on the reduced form models used to price credit derivatives. Nowak and Romaniuk (2013) apply TSIR models if the occurrence of the catastrophe does not depend on the financial market’s behavior. Zong-Gang and Chao-Qun (2013) derive a bond pricing formula in a stochastic interest rates environment with the losses following a compound nonhomogeneous Poisson process. Braun (2011) proposes a catastrophic swap pricing model representing the occurrence of catastrophes through a doubly stochastic Poisson process (Cox process) with a mean-reverting Ornstein-Uhlenbeck intensity. Lai, Parcollet and Lamond (2014) calculate the price of a catastrophe bond from a jump-diffusion process representing catastrophes, a three-dimensional stochastic process to represent the exchange rate, domestic and foreign interest rates, and the hedging cost for the currency risk.

Pérez-Fructuoso (2008) and Pérez-Fructuoso (2009) developed a new model for calculating the catastrophe loss index, whose value is the sum of the reported loss amount for each event. This variable is calculated as the difference between the total catastrophe’s severity and the key variable in the model named the incurred but not yet reported claims amount which is driven by a geometric Brownian motion with a constant claims reporting rate.

To develop a more precise expression of the catastrophic loss ratio, Pérez-Fructuoso (2016) develops an alternative model whose central hypothesis is that the incurred but not yet reported claims amount decreases proportionally to an exponential function, called the asymptotic claims reporting rate. The dynamics of this decrease is represented by a geometric Brownian motion.

Finally, and with the same objective as in the previous case, Pérez-Fructuoso (2017), models the decreasing linear dynamics of incurred but not yet reported claims amount, by means of an additive Brownian process or Ornstein-Uhlenbeck process.

Finally, Pérez-Fructuoso (2022) makes a comparison of the three models mentioned above. With the available data, it can be concluded that Ornstein-Uhlenbeck model is the one that fits better the real-life claims reporting process. However, we have also seen that the asymptotic model fits well the first two weeks after the catastrophes occurred.

A catastrophe loss index can be defined as the quotient between the total losses from catastrophes occurred over the period [0,T] (risk period) and a constant value p, whose definition depends upon the kind of index to be employed. The index value at maturity T′, LI(T′),with T′≥T, is defined as:

LI(T′)=L(T′)p (1)

where L(T′) is the total claims reported in T′ for all catastrophes occurring during the risk period.

L(T) is a random variable that depends on the following risk factors:

Under these assumptions, the numerator of the catastrophe loss index can be calculated as:

L(T′)=∑j=1N(T)Sj(T′) (2)

To obtain the value of (1), we can assume some hypotheses for our modelling.

It is well known in probability theory that; a Poisson process is a stochastic process that counts the number of events that happen in a certain period [0,T]. Therefore, we assume that N(T) is a Poisson process,

N(T)∼Poisson(λ)

where λ is the average number of catastrophic events per period [0,T]. The time between two different events in a Poisson process has an exponential distribution with the same parameter. Then, for our modelling:

tj+1−tj∼Exp(λ)

If Kj is the total amount of the catastrophe j and Sj(t) is for the reported loss amount at time t, we define Rj(t) as the incurred-but-not-yet-reported (IBNRL) loss amount at time t associated to the catastrophe j. Thus,

Sj(t)=Kj−Rj(t)

and the numerator of the catastrophe loss index (1) becomes:

L(T′)=∑j=1N(T)(Kj−Rj(T′)) (3)

In the following section, we develop the model that allows us to obtain an expression of Rj(t)to calculate the reported loss amount, and therefore the numerator of the catastrophe loss index triggering Cat Bonds. We consider the occurrence of a single generic catastrophe, K, occurring at moment τ∈[0,T], and we assume that its claims claim reporting process is S(t), and R(t) is its occurred but not yet reported loss process.

We define catastrophe severity as the sum of two random variables,

K=S(t)+R(t) (4)

were R(t) represents the incurred-but-not-yet-reported loss (IBNRL) amount at time t associated to the catastrophe occurred in τ and S(t) stands for the reported loss amount at time t related to a single catastrophe occurred in τ.

Once a catastrophe of severity K has occurred in τ, the claim reporting process associated with this single catastrophe is initiated lasting until the end of the maturity period. Empirical evidence shows that the intensity of claims reporting would seem to be greatest just after the occurrence of the catastrophe and decreases with time.

Then, the instantaneous claim reporting process is represented by a differential equation that describes the increase of reported loss amount proportional to the incurred-but-not-yet-reported loss amount, the main variable of the formal model,

dS(t)=α(t−τ)•R(t)dt (5)

were α(t−τ) is a real-value function named claim reporting rate function.

By differentiating equation (5) results,

dS(t)=−dR(t) (6)

and by substituting (6) into equation (5), the differential equation that describes the evolution of R(t) is obtained:

dR(t)=−α(t−τ)•R(t)dt (7)

Equation (7) shows that the incurred-but-not-yet-reported loss amount in t decreases with time according to the claim reporting rate function.

In order to capture the irregular behavior of the claim reporting process, we introduce a Wiener process (Brownian motion) into equation (7). This irregularity in the claim reporting process depends on the IBNRL amount. While there is still much to declare, the irregularity of the declarations is high. However, it decreases as the IBNRL amount does it too. To reflect this behavior, we will add a Wiener process with intensity σR(t), which is called a geometric Brownian motion. The result of adding it into our model (7) is the following stochastic differential equation,

dR(t)=−α(t−τ)•R(t)dt+σ•R(t)•dW(t) (8)

where α(t−τ) is the claim reporting rate function that represents the process’ drift, σ is a constant value that represents the process’ volatility and Wt is a standard Wiener process associated to the catastrophe.

A necessary condition for our modelling lies in the fact that σ must be a low value. Otherwise, if σ is large enough, the claim reporting rate might become positive and then the IBNRL amount will also grow, unlike our initial assumption.

Applying Itô’s Lemma (Friedman (1975); Malliaris and Brock (1991); Arnold (1974)) in equation (8), we obtain the expression for the incurred-but-not-yet-reported loss amount, R(t):

R(t)=K•exp(−∫0t−τα(s)ds−σ22•(t−τ)+σ•W(t−τ)) (9)

with the following boundary conditions:

From the relation defined between R(t) and S(t)defined in equation (4), we obtain S(t) as follows:

S(t)=K•[1−exp(−∫0t−τα(s)ds−σ22•(t−τ)+σ•W(t−τ))] (10)

Symmetrically to the incurred-but-not-yet-reported loss amount, the boundary conditions for the reported loss amount are:

Notice that is straightforward to see that if σ=0 we draw as a result the expression for both the incurred-but-not-yet-reported loss amount and the reported loss amount in a deterministic model:

R(t)=K•exp(−∫0t−τα(s)ds) and S(t)=K•[1−exp(−∫0t−τα(s)ds)] (11)

When the claim reporting rate is defined in a hybrid form, we assume that this rate is increasing up to a certain point in time, which we symbolize as tm, and thereafter this rate remains constant at the level α until the claims reporting process ends:

t_{m} \\ \end{matrix} \right.\ ]]> α(s)={αmtm•s0≤s≤tmαms>tm (12)

In proposed model to represent the claim process, the claim reporting rate is the variable that should be calculated for each set of real data. The model of claim reporting rate is defined through three parameters (αm,σ,tm). Thus, the global optimization procedure must simultaneously adjust all of them. The main goal of the global optimization problem is summarized in the following definition (Törn, 1991):

Given a function: f:M⊆ℜn→ℜ,M≠⌀,

for x∈M the value - \infty\ ]]> f*≔f(x*)>−∞is a global minimun, iff: ∀x∈M:f(x*)≤f(x)

Then x* is a global minimum point, f is called objective function, and the set M is called the feasible region. In this case, the global optimization problem has a unique restriction: the claim reporting rate must be positive. This restriction is included in the codification and all individuals are processed to become feasible ones. Then, despite this restriction, the solutions space does not have infeasible regions.

Evolutionary Algorithms are the term used to describe a broad set of algorithms that draw inspiration from the biological process of natural selection. Examples of evolutionary algorithms include genetic algorithms, genetic programming, evolutionary strategies, and differential evolution.  A notable feature that explains their success when applied to optimization problems is their ability to strike an appropriate balance between exploration and exploitation during the search process. Evolutionary Algorithms combine characteristics of both classifications of classical optimization techniques, volume-oriented and path-oriented methods. Volume-oriented methods (Monte-Carlo strategies, clusters algorithms) apply the searching process scanning the feasible region while path-oriented methods (pattern search, gradient descent algorithms) follow a path in the feasible region. A definition of a restricted search space of finite volume and the starting point is required to volume-oriented and path-oriented methods respectively.

Evolution strategies (ES) developed by Rechenberg (1971) and Schwefel (1981), have been traditionally used for optimization problems with real-valued vector representations. As Genetic Algorithms (GA) this are heuristics search techniques based on the building block hypothesis. Unlike GA, however, the search is basically focused on the gene mutation. This is an adapting mutation based on the likely the individual represents the problem solution. The recombination also plays an important role in the search, mainly in adapting mutation. ES are techniques widely used (and more appropriated than GA) in real-values optimization problems. ES offer practical advantages facing difficult optimization problems (Fogel, 1997). These advantages are: its conceptual simplicity, broad applicability, potentiality to use knowledge and hybridize with other methods, implicit parallelism, robustness to dynamic changes, capability for self-optimization and capability to solve problems that have no known solutions. A general ES is defined as an 8-tuple (Bäck, 1996):

ES=(I,Φ,Ω,Ψ,s,ι,μ,λ)

where:

In this work, the definition of the individual has been simplified: the rotation angles nα have not been considered, nα=0.

The mutation operator generates new individuals as follows:

σi′=σi•exp(τ′•N(0,1)+τ•N(0,1))

x⃗′=x⃗+σi′•N⃗(0⃗,1)

ES has several formulations, but the most common form is (μ,λ)-ES, where \mu \geq 1]]> λ>μ≥1, (μ,λ) means that μ-parents generate λ-offspring through recombination and mutation in each generation. The best μ offspring are selected deterministically from the λ offspring and replace the current parents. Elitism and stochastic selection are not used. ES considers that strategy parameters, which roughly define the size of mutations, are controlled by a “self-adaptive” property of their own. An extension of the selection scheme is the use of elitism; this formulation is called (μ+λ)-ES. In each generation, the best μ-offspring of the set μ-parents and λ-offspring replace current parents. Thus, the best solutions are maintained through generation. The computational cost of (μ,λ)-ES and (μ+λ)-ES formulation is the same.

In the case where the claim reporting rate is defined in a hybrid form, we calculate the integral of equation (9) in two different situations:

R(t) follows a log-normal distribution, where normal distribution parameters associated are:

lnR(t)∼N(lnK−(αm•(t−τ)2tm+σ22)•(t−τ);σ2•(t−τ))

lnR(t)∼N(lnK−(αm+σ22)•(t−τ)+αm•tm2;σ2•(t−τ))

As noted in section 3.1, the catastrophe loss ratio is defined by equation (1). For the case of catastrophe bonds, the indices used as triggers for indemnity payments are based on the accumulated losses to maturity associated with a single catastrophe. Therefore, we define the Bernoulli variable or indicator variable, δ, which is 0 if the catastrophe covered in the issue does not occur or 1 otherwise and the loss ratio is rewritten as:

LI(T′)=δS(T′)p={0siδ=0S(T′)p=K−R(T′)psi δ=1 (19)

When the claim reporting rate is defined through a hybrid model, the catastrophe loss index at maturity is given by the following expression:

LI(T′)=LI(tm)+LI(T′−tm)=

=δ•Kp•[(1−exp(−(αm•(tm−τ)2tm+σ22)•(tm−τ)+σ•W(tm−τ)))++(1−exp(e−(αm+σ22)•(T′−tm)+σ•W(T′−tm)•eαm•tm2))] (20)

In this work, we have a mathematical model that must be fitted to real data samples obtained from real catastrophes. Then, the problem could be defined as: "search for parameters that define a function that minimizes the error for each real data sample".

The proposed model follows equation (15) and (18), in both cases three parameters ( αm,σand tm)should be calculated to obtain the minimum distance to the real data distribution. The fitness function (evaluation of everyone over each series) is clearly the sum of the squared errors over the real data set but, in this case, due to the stochastic nature of the function, each individual must be evaluated by calculating the average value of the error over a large number of experiments repetitions. In this work, 10,000 evaluations have been performed to eliminate the randomness introduced by the Wiener process.

The type of recombination used in this work is the discrete recombination and the strategy ( μ+λ)-ES was used to select the individual to the next generation.

The algorithms used in this work have been developed by the authors in the R language and with the support of the cmaesr package (Bossek, 2021). The execution has been carried out on several computers Intel Core i5 4.8GHz with Windows 11 operating system. We have three computers with these characteristics to carry out the work. We performed 10,000 evaluations in each city to have the data available in one week. This number of evaluations is more than enough to ensure that the average value obtained is significant.

So, considering that we have data associated with seven cities and that each city has an execution time of 10000 evaluations of 1 week, the total execution time has been 3 weeks (first week three cities, second week three cities and third week the remaining city).

The parameters of the ES are summarized in Table 1. Besides, different runs were achieved changing the random seed1.

The real data to adjust model is related with floods occurred in several Spanish regions: Alcira (10/01/1991), San Sebastián (06/23/1992), Barcelona 1 (09/14/1999), Barcelona 2 (10/20/2000), Murcia (10/20/2000), Valencia (10/20/2000) and Zaragoza (10/20/2000). The data are temporal series of week incurred-but-not-yet-reported loss. For each disaster, we apply the optimization procedure for the model.

These data have been elaborated by the Technical and Reinsurance Department of the Consorcio de Compensación de Seguros (a public institution dependent on the Spanish Ministry of Economic Affairs and Digital Transition) to be applied exclusively in this research. The way in which they are presented, as a percentage of the weekly reported loss amount, means that they are not affected by the passage of time. That is, the data expressed as a percentage allows avoiding the time gap for its use at different moments in time. The amount of the catastrophe could be higher today, but the percentage declared in the first week would remain approximately the same (mainly because the population, public infrastructure and housing have not been modified in the affected areas after reconstruction).

Table 3 shows global results over the set of catastrophes for proposed model with evolutionary strategies, the accumulative quadratic error for each catastrophe and for the model. For each catastrophe, the model has different parameters (the best parameters for this real data distribution). Finally, the accumulative, the average and the standard deviation of the error for all catastrophes are shown.

For each event, the optimization procedure described above is applied. Table 4 shows the parameters of the proposed model associated with each of the catastrophes considered: