Catastrophe bond pricing for the two-factor Vasicek interest rate model with automatized fuzzy decision making

Catastrophe bonds are financial instruments, which enable to transfer the natural catastrophe risk to financial markets. This paper is a continuation of our earlier research concerning catastrophe bond pricing. We assume the absence of arbitrage and neutral attitude of investors toward catastrophe risk. The interest rate behavior is described by the two-factor Vasicek model. To illustrate and analyze obtained results, we conduct Monte Carlo simulations, using parameters fitted for real data on natural catastrophes. Besides the crisp cat bond pricing formulas, we obtain their fuzzy counterparts, taking into account the uncertainty on the market. Moreover, we propose an automated approach for decision making in fuzzy environment with relevant examples presenting this method.


Introduction
Nowadays, natural catastrophes occur more frequently than before. Additionally, they threaten densely populated areas Communicated (2011)]. This results in high values of damages. One can mention Hurricane Andrew (1992) with losses estimated at USD 30 billion [see, e.g., Muermann (2008)]. On the other hand, devastating floods are a significant problem in Europe.
Natural catastrophes losses have a negative impact on financial stability of insurers, e.g., after Hurricane Andrew more than 60 insurance companies became insolvent [see, e.g., Muermann (2008)]. This is caused not only by the huge value of losses, but also the classical insurance mechanisms [see, e.g., Borch (1974)] are not suitable for losses caused by natural disasters. The classical insurance approach applied to the construction of the insurer's portfolio can lead directly to bankruptcy of this enterprise [see, e.g., Ermoliev et al. (2001); Romaniuk and Ermolieva (2005)].
One of the methods of solving the mentioned above problems is application of other financial and insurance mechanisms. Such instruments are, among others, catastrophe bonds [also called cat bonds, see, e.g., George (1999); Nowak and Romaniuk (2009)]. They transfer risk connected with natural catastrophe losses (i.e., the insurance risk) to financial markets.
The number of papers dedicated to catastrophe bonds pricing is relatively low. A significant problem is the incompleteness of the financial market on which catastrophe instruments are traded. In many papers, simplified models or behavioral finance methods are used. In turn, many stochastic approaches do not take into account the change of probability measure in the martingale method or apply the utility function, which choice can be an additional problem in practice. An interesting stochastic approach proposed by Vaugirard in Vaugirard (2003) overcame the problem of the incompleteness of the financial market. For more detailed references we refer the reader to Nowak and Romaniuk (2013b).
In this paper, we solve the catastrophe bond pricing problem under the assumption of absence of arbitrage on the financial market. We derive and prove the cat bond valuation expressions in crisp case. Then, assuming that some parameters of the model cannot be precisely described, we obtain fuzzy counterparts of the crisp pricing formulas and we apply them to an automatized method of decision making. We illustrate our theoretical results by numerical examples. Our approach in this paper can be treated as an essential extension of Vaugirard's method. In analogous catastrophe bond valuation formulas considered earlier [see Nowak et al. (2012); Nowak andRomaniuk (2010b, 2013b)], one-factor models of the spot interest rate r t were used. A new model of r t proposed in this paper is the two-factor Vasicek model (introduced by Hull and White) with a stochastic process describing deviation of the current view on the long-term level of r t from its average view. Additionally, in comparison with the mentioned approaches, we introduce more detailed formulas describing prices of cat bonds (see Theorem 2) for arbitrary time moment t, on the basis of Dynkin's Lemma. For this purpose, we apply zero-coupon bond pricing formulas from Munk (2011). However, the proof of Theorem 2 required considering the case of a r = a ε (see Lemma 1, case b), not considered in Munk (2011). Furthermore, contrary to Nowak et al. (2012), Nowak and Romaniuk (2010b), we conduct Monte Carlo simulations, applying real parameters for r t and distributions of catastrophe losses. We also apply a non-homogeneous Poisson process for quantity of catastrophe events. As mentioned previously, apart from the standard valuation formulas, we obtain their fuzzy counterparts, taking into account the uncertainty on the market. We propose an automated approach for decision making and present examples to illustrate it for a given set of market parameters.
As noted above, applying fuzzy set theory, we consider different sources of uncertainty, not only the stochastic one. In particular, we take into account that volatility parameters and the correlation coefficient, used for description of r t , are determined by fluctuating financial market and often the uncertainty does not have a stochastic character. Therefore, we apply fuzzy counterparts of some market parameters. As result, price obtained by us has the form of a fuzzy number, which can be used for investment decision making. Similar approach was applied in the case of options in Wu (2004) and Nowak and Romaniuk (2010a, 2013c, 2014b, where the Jacod-Grigelionis characteristic triplet [see, e.g., Nowak (2002); Shiryaev (1999)] was additionally used.
The literature devoted to soft computing methods is very rich. Books (Vasant et al. 2012;Vasant 2013aVasant , 2014Vasant , 2015 provide a wide overview of soft computing algorithms and their applications. We mention only certain selected examples of the presented approaches, focusing on financial and economic problems. Various optimization methods using soft computing such as fuzzy logic, neural networks, genetic algo-rithms, and the theory of chaos, in business, finance and economics are discussed in Dostál (2012). An interesting example presented in Dostál (2012) is application of theory of chaos to simulation of market price of a share on the stock market. In Beynon and Clatworthy (2013), the problem of understanding the relationship between company stock returns and earnings components is discussed. The authors apply the classification and ranking belief simplex (CaRBS) and a development of this data analysis technique, called RCaRBS. In turn, the author of Gordini (2014) uses a genetic algorithms approach for small enterprise default prediction modeling. A widely applicable multi-variate decision support model for market trend analysis and prediction, based on a time series transformation method in combination with Bayesian logic and Bayesian network, is presented in Mršić (2014). In Besbes et al. (2015), the authors propose a twophase mathematical approach, involving the application of a stochastic programming model, to effective product-driven supply chain design. In Vasant (2013b), a new and interesting approach to the industrial production planning problems is discussed. Finally, various generalized soft methods are applied in many other scientific areas, see, e.g., Kowalski et al. (2008), Charytanowicz (2005, 2010).
Section 2 contains lists of symbols, expressions and operators used in the paper. In Sect. 3, we discuss basic notions and mechanisms connected with catastrophe bonds. In Sect. 4, catastrophe bond pricing formulas are proved for the twofactor Vasicek interest rate model and a necessary theoretical base is presented. Section 5 is devoted to cat bond pricing in fuzzy environment. In Sect. 6, we propose an automated method of decision making with application of the fuzzy valuation expressions. Section 7 is dedicated to algorithms used in Monte Carlo simulations. These algorithms are then applied in Sect. 8 during analysis of some numerical examples of cat bond pricing and automatized decision making.

Catastrophe bonds
A single catastrophic event, e.g., an earthquake or a hurricane, could result in damages worth of billions of dollars. Therefore, such event could cause problems with reserves for many insurers or even bankruptcy of these enterprises [see Cummins et al. (2002)]. However, daily fluctuations on worldwide financial markets reach also tens of billions of dollars. Then, new instruments were introduced, which could transfer risk from insurance markets to financial markets [see, e.g., Nowak (1999)].
A catastrophe bond (Act-of-God bond, cat bond) is an example of such instruments aimed at "packaging" risks into a form of tradable assets [see Cox et al. (2000); George (1999); Nowak and Romaniuk (2009)]. The payment function of the cat bond depends on additional random variable called triggering point, which is connected with occurrence or other properties of specified type of natural catastrophe (like the issuer's actual losses-e.g., losses from flood, insurance industry index, real parameters of catastrophe-e.g., magnitude of earthquake, etc.). If triggering point occurs, then the structure of payments is changed. Other parameters like region and time interval for catastrophic event are described in detail for catastrophe bond. The payments for cat bonds usually depend also on interest rates.
For example, the A-1 bond issued by USAA in 1997 was connected with losses caused by hurricane on the east coast of USA between July 15, 1997 and December 31, 1997. If the value of losses had been more than $ 1 billion, the coupon of the bond would have been lost.

Cat bond pricing
This section is devoted to catastrophe bonds pricing in the crisp case, i.e., valuing the catastrophe instruments under the assumption that all the model parameters are crisp. At the beginning, we present selected elements of stochastic analysis, describe the stochastic model of catastrophe losses and introduce Brownian motion, which is used for description of the risk-free spot interest rate. We also define two types of catastrophe bonds and specify assumptions concerning the considered financial market. Subsequently, we discuss the two-factor Vasicek model of the spot interest rate and complete the zero-coupon bond valuation expression for this interest rate model. Finally, we introduce and prove the catastrophe bonds valuation formulas.
We introduce necessary notations and definitions. Let ( , F, P) be a probability space. We assume that time horizon has the form 0, T , where T < ∞, and we fix a moment T ∈ 0, T .
be a sequence of independent random variables with the same distribution and finite second moment. We treat U i as the value of losses during ith catastrophic event. To describe the aggregated catastrophe losses till moment t, we will use the compound Poisson process given by formulã where Moreover, we assume that F = F T , Let (B t ) t∈[0,T ] denote the banking account, satisfying the standard equation where r is the risk-free spot interest rate. We will denote by the symbol B (t, T ) the price at the time t of a zero-coupon bond with the maturity date T ≤ T and with the face value equal to 1 .
We attempt to introduce definitions of two types of catastrophe bonds with different payoff structures. Let 0 < K 1 < · · · < K n , n ≥ 1, be levels of catastrophic losses and let The constants w i will be called payment decrease coefficients.
In this section, we denote by ∧ the two-argument minimum function.
Let τ i : → 0, T , 1 ≤ i ≤ n be a sequence of stopping times defined by the formula Definition 1 By the symbol IB s (T, Fv), we denote catastrophe bond with the face value Fv, the date of maturity and payoff T , satisfying the following conditions: (a) If the catastrophe does not occur in the period [0, T ] (i.e., τ 1 > T ), the bondholder is paid the face value Fv.
Clearly, the catastrophe bond IB s (T, Fv) payoff function ν IB s (T,Fv) can be written in the form: where I is the indicator function.
In the reminder of this section, we will use the following function of the variable T : where i are cumulative distribution functions of τ i .
To define the second type of cat bond, we introduce an additional constant K 0 , such that 0 ≤ K 0 < K 1 .
Definition 2 By the symbol IB p (T, Fv), we denote catastrophe bond with the face value Fv, the date of maturity and payoff T and the payoff function of the form The payoff function of the cat bond IB p (T, Fv) is a piecewise linear function of lossesÑ T . If the catastrophe causing significant level of losses does not occur (i.e.,Ñ T < K 0 ), the bondholder receives the payoff equal to its face value Fv. In case when the aggregated losses are not less then K n , the bondholder receives the payoff equal to Fv arbitrage-free family of zero-coupon bond prices with respect to r , if the following conditions hold: There exists a probability measure Q, equivalent to P, such that for each T ∈ 0, T , the process of the discounted price of the zero-coupon bond is a martingale with respect to Q.
Letλ t = λ 1t ,λ 2t , . . . ,λ nt be the n-dimensional market price of risk process. In our model, we assume thatλ 1t = λ 1 ,λ 2t = λ 2 , . . . ,λ nt = λ n are constants. The following Radon-Nikodym derivative defines a probability measure Q equivalent to P: where . is the Euclidean norm in R n . The family of zero-coupon bond prices B(t, T ), where 0 ≤ t ≤ T ≤ T is arbitrage free with respect to r for Q defined by (2).
In what follows we assume that investors are neutral toward the natural catastrophe risk. Additionally, we introduce the following notations: From the independence of the payoff function from W [compare Nowak et al. (2012) and Nowak and Romaniuk (2010b)], it follows that

The two-factor Vasicek model
We now introduce the pricing formulas for the catastrophe bonds under the assumption of the two-factor Vasicek model of the spot interest rate. In the mentioned model, the interest rate is described by the system of stochastic differential equations, which in the form proposed by Hull and White [compare Hull and White (1994)] is given by The process (ε t ) t∈ [0,T ] represents the deviation of the current view on the long-term level of the interest rate (r t ) t∈ [0,T ] from the average view. The parameter ρ ∈ [−1, 1] is the correlation coefficient between changes in the interest rate and changes in the process ε. We assume that all the parameters (except the correlation coefficient ρ) in the above equations are positive. The proposed model belongs to the class of Gaussian models.

Cat bonds valuations formulas
In catastrophe bond pricing formulas, we will apply the zero-coupon bond valuation expression from Munk (2011). However, since the author of the mentioned book did not take into account the case of a r = a ε , we will use the following Theorem 1 (Theorem 8.1 adapted to the two-factor affine interest rate models) from Munk (2011) to complete the zerocoupon bond pricing expression. The complete expression of B(t, T ) will be presented in Lemma 1.
Let us assume that the spot interest rate process r under Q has the form for the dynamics of the vector of state variables given by t andW 2 t are independent Q-Brownian motions. We additionally assume that for all possible values of the vector of state variables.

Theorem 1
In an admissible affine model, the price of a zero-coupon bond has the exponential-affine form where the deterministic functions a, b 1 , b 2 satisfy the system of ordinary differential equations Functions a, b 1 , b 2 in the above theorem are auxiliary and they are used for description of the price of a zero-coupon bond.
In the following lemma, we assume that for a constant λ (compare the proof of Theorem 2). To shorten formula describing B (t, T ), we introduce additional auxiliary functions: B a r , B a ε , B a r +a ε , C, D, E, G.
Lemma 1 Let the risk-free spot interest rate be described by the two-factor Vasicek model. Then, the zero-coupon bond pricing formula has the following affine-exponential form: where for τ ∈ [0, T ]: and and Proof For the considered two-factor Vasicek model, the system of the ordinary differential equations formulated in Theorem 1 has the form with the initial conditions a (0 The solution of (18) in the case of a r = a ε has the form (7). Moreover, the solution of (19) with b 1 and b 2 of the form (6) and (7), respectively, is given by (5). The zero-coupon bond pricing formula obtained this way was presented and proved in Munk (2011).
In the case of equal coefficients a r , a ε , the solution of (18) has the form (12). We substitute formulas (11) and (12) into equation (19). Then, we use the equality to obtain the solution of (19). Applying standard integral operations for right side of (19), we receive the form of a (τ ) described by (10) with functions C, D, E, G given by (13)-(16) and, finally, the zero-coupon bond pricing formula of the form (4).
Theorem 2 Let IB s (t) and IB p (t) be prices of the catastrophe bonds IB s (T, Fv) and IB p (T, Fv) at moment t ∈ [0, T ] for the risk-free spot interest rate described by the two-factor Vasicek model. Then and where functions a, b 1 , b 2 are described by formulas (5)- (7) in the case of a r = a ε and by (10)-(12) for a r = a ε .
Proof We assume that the market price of risk connected with the deviation of the current view on the long-term level of r t is equal to zero. Then, if λ 1 = λ is a fixed constant, λ 2 satisfies the following equation: After change of measure according to the formula (2), the Eq. (3) has the form: and, therefore, from Lemma 1, it follows that where a, b 1 , b 2 are described by formulas (5)-(7) in the case of a r = a ε and by (10)-(12) for a r = a ε . Moreover, for each t ∈ [0, T ] Fv) . The above formula follows from Dynkin's Lemma applied to the π -system where c t = exp − T t r u du . For more details, we refer the reader to the one-factor affine interest rate case considered in Nowak and Romaniuk (2013a). Applying the equality (22) to formula (23), we obtain (20) and (21).

Fuzzy approach to catastrophe bonds pricing
In this section, we present and prove the catastrophe bonds pricing formulas, assuming that some model parameters are not precisely known and they are described by fuzzy numbers. We recall basic elements of fuzzy number theory and fuzzy arithmetic. Using Zadeh's extension principle, we prove fuzzy counterparts of catastrophe bonds valuation expressions introduced in the previous section. Additionally, we derive analytical formulas for α-level sets of the cat bonds prices. Finally, we discuss a method of computation of values of their membership functions.

Basic notions and facts
Now we recall some basic notions and facts concerning fuzzy numbers and fuzzy and interval arithmetic.
LetÃ be a fuzzy subset of the set of real numbers R. We denote by μÃ its membership function μÃ : For a fuzzy numberã [see, e.g., Nowak and Romaniuk (2014b) Let be a binary operator ⊕, , ⊗ or between fuzzy numbersã andb, corresponding to the related binary operator • between real numbers: +, −, × or / via the extension principle. Let int be a binary operator ⊕ int , int , ⊗ int or int between two closed intervals [a, b] and [c, d], defined by where • is the corresponding to int operator between real numbers +, −, × or / (if the interval [c, d] does not contain zero in the last case).
Letã,b be fuzzy numbers. Then,ã b is also a fuzzy number and the following equalities hold: if α-level setb α does not contain zero for all α ∈ [0, 1] in the case of . A triangular fuzzy numberã = (a 1 , a 2 , a 3 ) is a fuzzy number with the membership function of the form We denote by F (R) the set of all fuzzy numbers. In the further part of this section, we will use the following proposition, proved in Wu (2004).
We will also use fuzzy random variables. For their definition, we refer the reader to Puri and Ralescu (1986).

Pricing formulas
In practice, some parameters of the financial market (e.g., the market price of risk, the correlation coefficient or the volatility parameters) are not precisely known. Many authors noticed that the uncertainty of the market parameters does not have a stochastic character [see, e.g., Wu (2004)]. Since the parameters are determined by the market which fluctuates from time to time, it is unreasonable to choose their fixed values, obtained from historical data, for later use in pricing formulas. Therefore, to estimate values of the mentioned parameters one can use knowledge of experts, regarding the market price of risk, the correlation coefficient, and the volatility parameters as fuzzy numbers. One can ask experts for forecasts of the parameters values, transfer them into triangular fuzzy numbers, and then use them for estimation of the parameter.
In the remainder of this section, we assume that the abovementioned parameters are fuzzy numbers of the triangular form. In particular, we introduce fuzzy numbersσ r ,σ ε ,λ 1 , ρ and fuzzy random variablesr t ,ε t , t ∈ [0, T ] , in place of their real counterparts σ r , σ ε , λ 1 , ρ and r t , ε t , t ∈ [0, T ]. We assume thatσ r andσ ε are non-negative fuzzy numbers, i.e., their membership functions are equal to 0 for all negative arguments.
We denote by the set of symbols = {L , U } and introduce the operator : → by: L = U and U = L.
In the following theorem, we will prove the fuzzy valuation formulas for the two-factor Vasicek interest rate model. Simultaneously, we will correct a simplified version of these formulas from Nowak and Romaniuk (2014a), where only fuzzy volatility parameters were assumed.

Computational methods
Apart from the forms of α-level sets of the catastrophe bond price, it is necessary to compute the value of its membership function μĨ B(t) (c) for arbitrary c. We fix a time moment describes the membership function ofĨ B (t).
We will use the above formula for computation of the membership function μĨ B(t) , applying the method proposed for pricing European options in Wu (2004). From the pricing formulas (see Theorem 3), it follows that the functions Problem (OP1) can be rewritten to the form: To solve (OP2), it suffices to consider the following three cases: (1), then one can solve the following relaxed optimization problem: (1), then one can solve the following relaxed optimization problem: Problems (OP3) and (OP4) can be solved using bisection search [see Wu (2004)].
In the further part of the paper, we will use the following functions ofĉ: One can check that

Automatized method of decision making
The obtained catastrophe bond pricing formulas can be used for investment decision making. We apply an approach similar to the one used in Piasecki (2014) in another context. We obtain a fuzzy set of possible decisions. A decision maker can choose one of the decisions from its α-level set for a sufficiently high value of α. Let t ∈ [0, T ]. For a given α (e.g., α = 0.9), the α-level set ofĨ B (t) can be treated by a financial analyst as the interval of the cat bond prices. The financial analyst can choose any value from this interval as the catastrophe bond price with an acceptable membership degree. Such an interval can be a very useful tool for investment decision making.
In particular, apart from a straightforward analysis of αlevel sets ofĨ B (t), the financial analyst can consider the following set of possible decisions: The advice choice function : R 2 → 2 V is defined as follows: where the symbol AND denotes logical conjunction. Then, by applying the Zadeh extension principle, we obtain the extended advice choice functioñ The membership functionl of˜ Ĩ B (t) ,Ĉ t is defined by formulas: The financial analyst can choose one of the decisions from the α-level set˜ Ĩ B (t) ,Ĉ t α for a sufficiently high value of α.

Algorithms used in simulations
We start from general remarks concerning algorithms which are used in catastrophe bond pricing (see Sect. 7.1 for the crisp case and Sect. 7.2 for the fuzzy case) and in automatized method of decision making (see Sect. 7.3). Then, these algorithms are applied to obtain numerical results in some exemplary cases considered in Sect. 8.

Cat bond pricing-crisp case
To find the price of the cat bond in the crisp case, we apply Monte Carlo simulations to the results obtained in Theorem 2.
From (20) and (21), it is seen that numerical complexity of the considered model and the related pricing formula depends on two factors: the expected value of the payoff function e t or f t and the discounting factor related to adopted in this paper two-factor Vasicek model. In the following, Monte Carlo approach is applied to approximate e t or f t . The error of such estimation measured by standard deviation is proportional to 1/ √ n S , where n S is the number of simulations. But during analysis of the computational complexity, the form (1) of the processÑ t should be also taken into account, because the numerical efficiency of simulation of the single trajectory ofÑ t depends on both types of the Poisson process N t and the random distribution of U i . We assume that non-homogeneous Poisson process (NHPP) is applied for modeling quantity of losses and that the value of single loss is modeled by lognormal distribution or Pareto distribution. These distributions are commonly used in modeling of risk events in insurance [see, e.g., Chernobai et al. (2006)]. The parameters of the applied distributions were fitted in Chernobai et al. (2006) for real data describing natural catastrophic events in the United States provided by the ISO's (Insurance Service Office Inc.) Property Claim Services (PCS). The sinusoidal intensity function for NHPP given by the formula results in a better calibration than simpler homogeneous Poisson process as discussed in Chernobai et al. (2006). Because of the cyclic form (44) of κ(t), thinning method [see, e.g., Law (2007)] is applied during simulations. Rejection step in this approach results in some loss of the samples generated from the relevant homogeneous Poisson process. For the most straightforward upper bound of (44) given by a + b2π , the average acceptance rate for the interval [0, T ] is equal to Algorithms for generation of variables from lognormal or Pareto distribution are numerically very fast, like e.g., Marsaglia-Tsang ziggurat algorithm for embedded normal distribution [see Marsaglia and Tsang (2000) for details]. In the case of Pareto distribution, the direct method of inversion of cdf is used. Therefore, to obtain one random variable from the target distribution, only one output from PRNG is sufficient.
Of course other types of NHPP than the formula (44) or other loss distributions could be easily incorporated into the simulation approach. If more complex pdf is necessary, then even Markov Chain Monte Carlo method should be applied.
The whole procedure for finding estimator for e t or f t constitutes the following algorithm:

Algorithm 1
Input Parameters of: NHPP process, distribution of losses, model of payment; number of simulations n S .
Step 1 Sample trajectory of the processÑ t based on times of the losses generated according to κ(t) and values of the losses given by iid sample U 1 , U 2 , . . . generated from the distribution of losses.
Step 2 Calculate the relevant payoff function ν IBs (T,F v) or ν IBp(T,F v) using the output from the step 1.
Step 3 Store the evaluated value of payoff.
Step 4 Return to the step 1, unless the necessary number of simulations n S is acquired.
Step 5 Calculate the average value of the payoff based on the values stored during the step 3. Output Monte Carlo estimator of the expected value of the payoff function e t or f t .
In this paper, we assume that the model of interest rates is described by the two-factor Vasicek model. In the following, the parameters for this model fitted in Puhle (2008) are used. Easily seen, two-factor model may be better calibrated to the real-life data than in simpler one-factor approach. Because in Theorem 2, the analytical form of the discount factor is given, then the relevant cat bond price could be directly evaluated without the necessity of subsequent simulations after using Algorithm 1. The price IB s (t) or IB p (t) is straightforward multiplication of the previously estimated expected value of payoffs and the relevant discount factor.
The whole procedure is very fast. We use Intel Core 2 Quad (2.5 GHz) with 4GB RAM, VS Studio 2013 Express compiler, GSL library for random number generation and Mersenne Twister as primary PRNG. Even without sophisticated optimization (like multithreading), finding cat bond price estimator based on sampling 1000000 trajectories takes less than 10 seconds.

Cat bond pricing-fuzzy case
Now we focus on the cat bond pricing in the fuzzy case. In the following, especially in Sect. 8.2, only triangular fuzzy numbers are considered, but the introduced numerical approach may be also used for other types of LR numbers.
As it is seen from (28) and (29), similarly to the crisp case, the complexity of the considered model and the relevant pricing formulas depends on the (crisp) expected value e t or f t , and the α-level set related to the discounting factor. The estimator of e t and f t could be directly found using the Algorithm 1, regardless of the considered value α.
To approximate the fuzzy cat bond price, we fix the value of α and then find the relevant interval using formulas (28) or (29) depending on the applied payment function. If α is gradually set to subsequent values, starting from some initial value α 0 ∈ [0, 1] up to upper bound α 1 (where α 1 ∈ [α 0 , 1]) with given increment α > 0, then the obtained intervals putting on one another form the final output-approximation of the fuzzy cat bond price. It leads us to the following approach: Algorithm 2 Input Parameters of: NHPP process, distribution of losses, model of payment, two-factor Vasicek model; number of simulations n S ; α-level values: α 0 , α 1 , Δα.
Step 1 Using Algorithm 1, estimate the expected value of the payoff function e t or f t .
Step 2 Set (crisp) values of the parameters a r , a ε , b r , r 0 , ε 0 .
Step 5 Based on the crisp parameters from the step 2, the α-level sets obtained in the step 4 and the value of the estimator found during the step 1, calculate the interval of the cat bond price [g(α), h(α)] using (28) or (29), respectively.
Step 6 Store the obtained data (i.e. the α-level set) in order to approximate the output (i.e. the fuzzy cat bond price).
Step 8 Return to the step 4, unless α > α 1 . Output The approximated fuzzy cat bond price formed by intervals stored during the step 6.
Apart from the Monte Carlo simulations necessary during the step 1, the computational complexity of the considered fuzzy model directly depends only on method of approximation of the final fuzzy cat bond price. The analytical formulas for endpoints of the α-level sets of prices are given by Theorem 3. Therefore, for n 0 possible values of α approximating the fuzzy output, only 2n 0 m endpoints for α-level sets of the parameters should be used (where m is number of fuzzy parameters in the model). The above procedure is very fast. Without the first step, using the same platform as in Sect. 7.1, approximation of the fuzzy price for n 0 = 100 takes 0.054 s.
It should be noted that availability of the fuzzy analytical formulas provided by Theorem 3 is important advantage of the pricing method proposed in this paper. Otherwise, set of even n 0 2 m possible endpoints for α-level sets of the parameters should be considered to find fuzzy output.

Evaluation of the membership values
In this section, we consider the approaches to evaluate membership value α for the current market price of the cat bondĈ t which is also useful for automatized decision making introduced in Sect. 6.
As explained in Sect. 5.3, special optimization procedure based on bisection search is used to evaluate membership value α for the given current market price of the cat bondĈ t , i.e., is used to find μĨ B(t) Ĉ t solving the problem (OP2) [see also Wu (2004) for detailed description of the relevant algorithm]. In this approach, only the currently used α varies from step to step of the algorithm, so Monte Carlo estimator of e t or f t has to be found via Algorithm 1 only once at the beginning of the whole procedure. For the same platform as in Sect. 7.1, evaluation of the membership value for the givenĈ t with tolerance = 0.00000001 takes less than 0.01 seconds.
The evaluated membership value could be also used as a source of investment decisions, i.e., to obtain membership functions for the extended advice choice function˜ introduced in Sect. 6.
To calculate the membership functions for various possible decisions from the set V for the fixed current market priceĈ t , it is sufficient to evaluate values βĨ B(t) Ĉ t and δĨ B(t) Ĉ t used in the formulas (39)-(43). Then, we have the following procedure:

Algorithm 3
Input Parameters of: NHPP process, distribution of losses, model of payment, two-factor Vasicek model; number of simulations n S ; tolerance ; current market priceĈ t .
Step 1 Find g(1) and h(1) using Algorithm 2 for the fixed single value α 0 = 1. Store the estimators of e t or f t evaluated in the step 1 of Algorithm 2.
Output Values βĨ B(t) Ĉ t and δĨ B(t) Ĉ t for the fixed C t which are used to calculate the membership functionsl (B) ,l (A) ,l (H) ,l (R) ,l (S).
Easily seen, even for variousĈ t it is necessary to use Monte Carlo simulations to evaluate e t or f t only once and then this estimator could be reused. Therefore, the complexity of the main part of Algorithm 3 (apart from the step 1) depends only on the method to solve the problem (OP2), i.e., the bisection search.
To find approximation of the whole membership function for the extended advice choice function˜ for various α or only for the single fixed value of α, the approach similar to Algorithm 3 could be adopted for the whole range of current prices. As mentioned previously, βĨ B(t) Ĉ t and δĨ B(t) Ĉ t should be found for [Ĉ g t ,Ĉ h t ] starting from some initial valuê C g t < g(0) and stopping withĈ h t > h(0) with fixed incrementation Ĉ t > 0. Then, we have the following procedure:

Algorithm 4
Input Parameters of: NHPP process, distribution of losses, model of payment, two-factor Vasicek model; number of simulations n S ; tolerance ; market priceŝ C g t ,Ĉ h t , ΔĈ t .
Step 1 Find g(1) and h(1) using Algorithm 2 for the fixed single value α 0 = 1. Store the estimators of e t or f t evaluated in the step 1 of Algorithm 2.
Step 3 Go to the step 2 of Algorithm 3.
Step 4 Evaluate membership functions for the recommendations based on the output of Algorithm 3.
Step 6 Return to the step 3, unlessĈ t >Ĉ h t . Output Values βĨ B(t) Ĉ t and δĨ B(t) Ĉ t used to calculate the membership functionsl for the set of the possible decisions V for the whole interval of prices [Ĉ g t ,Ĉ h t ].
Once again, the most numerically complex part (i.e., Monte Carlo simulations in the step 1) is invoked only once, regardless of the number of points evaluated in the interval Ĉ g t ,Ĉ h t . For each of this point, the bisection solution in Algorithm 3 is applied only once. Using the same platform as in Sect. 7.1 and for n 0 = 100, the whole procedure (apart from finding estimators for e t or f t ) takes 0.26 s.

Numerical examples and analysis
Based on algorithms introduced in Sect. 7, some numerical examples of cat bond pricing (in both the crisp and the fuzzy environments, see Sects. 8.1 and 8.2, respectively) and automatized decision taking (in the fuzzy case, see Section 8.3) are considered. In each experiment, we conduct n S = 1000000 simulations. The relative percentage change of the cat bond price ( IB p (0)) and the relative percentage change of ρ ( ρ) are compared to the reference cat bond price and ρ from Example I.1

Cat bond pricing-crisp case
Algorithm 1 from Sect. 7.1 with formulas from Theorem 2 is used to evaluate the cat bond price in the crisp case. We assume that the face value of the bond in each experiment is set to 1 (one monetary unit assumption) and that the payment function is described by Definition 2. Then, the (crisp) cat bond price IB p (0) is given by (21).

Example I.1
As mentioned previously, model of the process of the losses and model of the interest rates are fitted to the real-life data. Then, the parameters of the intensity function (44) are equal to a = 30.8750, b = 1.6840 and c = 0.3396 [see Chernobai et al. (2006)]. The parameters of the lognormal distribution are set to μ LN = 17.3570 and σ LN = 1.7643 [see Chernobai et al. (2006)]. The two-factor Vasicek model is described by the data for term structure of interest rates and a r = 0.2591, b r = 0.0205, σ r = 0.0073, a ε = 0.8274, σ ε = 0.0219, ρ = 0.6, r 0 = 0.025, ε 0 = 0 [see Puhle (2008)].
We set n = 2 and w 1 = 0.25, w 2 = 0.5. The relevant levels of losses (also called triggering points) are connected with surpassing the limits given by quantiles of the cumulated value of losses described by the NHPP process (number of losses) and lognormal distribution (value of each loss). Such xth quantile is denoted by Q LN−NHPP (x). Similar approach was considered in Nowak and Romaniuk (2013b). We assume that K 0 = Q LN−NHPP (0.5), K 1 = Q LN−NHPP (0.75) and K 2 = Q LN−NHPP (0.95). Then, IB p (0) is estimated as 0.830485.
In the following, this cat bond price is treated as the reference value for subsequent comparisons. In the next examples, we change values of one or more parameters of the model to verify if these alternations have significant impact on the calculated cat bond prices.

Example I.2
We analyze the relation between the cat bond price and the value ρ for two-factor Vasicek model, i.e., we find the prices for ρ ∈ [0.1, 0.9]. The other parameters are the same as in Example I.1. Then, only the correlation coefficient for interest rate model is changing. Therefore, we could analyze the situation, e.g., if this parameter is mismatched during the estimation or it will vary in near future compared to the historical data. The outcomes could be also compared with the fuzzification of ρ considered in Sect. 8.2, which reflects other sources of uncertainty. The exact cat bond prices IB p (0) are given in the third row of Table 1. In the fourth row, the relative percentage change of the actual cat bond price IB p (0) compared to the reference value from Example I.1 is calculated. The relative percentage change ρ of the actual ρ compared to the reference value ρ = 0.6 is given in the second row. As it is seen from Fig. 1, the price as the function of ρ is increasing and linear, but the variability of the prices is almost negligible ( IB p (0) is less  The relative percentage change of the cat bond price ( IB p (0)) and the relative percentage change of a ( a) are compared to the reference cat bond price and a from Example I.1 than 0.0007 % even for ρ = 0.1). It should be noted that such numerical analysis is useful and straightforward, because the function of the price depends on ρ in rather complex way.

Example I.3
We analyze the relation between the cat bond price and the value a for the intensity function (44), i.e., we find the prices for the wide interval a ∈ [23,40]. Then, the overall intensity of the number of losses is lesser or higher than in the case of a fitted to the historical data in Chernobai et al. (2006).
Once again, such difference may be caused, e.g., by estimation error or historical data irrelevant for the future. The cat bond price as the function of a is decreasing and almost linear (see Fig. 2 and the third row of Table 2 for some exact values). But in this case the relative percentage changes of price are noticeable (decrease up to almost 9 % for the relative growth of a which is equal to about 26 %). Because of the applied intensity function (44)  The relative percentage change of the cat bond price ( IB p (0)) and the relative percentage change of b ( b) are compared to the reference cat bond price and b from Example I.1 is expected, but the detailed nature of such phenomena has to be numerically determined.

Example I.4
We analyze the relation between the cat bond price and the value b in the intensity function (44), so we find the prices for the wide interval b ∈ [0.5, 2.6]. Then, the intensity of number of losses given by the cyclic part in (44) fluctuates less or more than in the case of historical data from Chernobai et al. (2006), so the number of catastrophes is almost on the same level constantly or more significantly varies from period to period. The cat bond price as the function of b is decreasing and almost linear (see Fig. 3; Table 3) and the relative percentage changes are noticeable, but also lower compared to the previous example (about 6 % decrease of the price if b is increased more than 50 % compared to the reference values from Example I.1). As in the case considered in Sect. 8.1.3, the exact behavior of this function in such setting could be directly analyzed using simulations.

Example I.5
We analyze the relation among the cat bond price and the volatilities σ r and σ ε in the interest rate model, so we find the prices for σ r ∈ [0.001, 0.015] and σ ε ∈ [0.005, 0.04]. The behavior of these parameters is also considered in the fuzzy case in Sect. 8.2. The other parameters are the same as in Example I.1. Therefore, the volatilities vary compared to the historical data from Puhle (2008).
The cuts of the graph (see Fig. 4) for both volatilities (i.e., if one of them is set as constant) seem to be hyperbolic functions. The estimated values of the cat bond prices for σ ε = 0.0219 and σ r = 0.0073 are given in Tables 4 and  5  The relative percentage change of the cat bond price ( IB p (0)) and the relative percentage change of σ r ( σ r ) are compared to the reference cat bond price and σ r from Example I.1 The relative percentage change of the cat bond price ( IB p (0)) and the relative percentage change of σ ε ( σ ε ) are compared to the reference cat bond price and σ ε from Example I.1 relative changes of prices are comparable in scale for both the volatilities.

Example I.6
We apply generalized Pareto distribution (GPD) instead of lognormal distribution for the value of the single loss. The parameters of the GPD are set to ξ GPD = 0.8090, β GPD = 5.340 × 10 7 (see Chernobai et al. (2006)), so they are estimated in the same manner and from the same real-life data as in the previous cases. To enable comparisons, we apply approach similar to Example I.1, so the two-factor Vasicek model is described by the same parameters and we also set n = 2 and w 1 = 0.25, w 2 = 0.5. But the relevant triggering points are given as quantiles of the cumulated value of losses described by the NHPP process and GPD, instead of lognormal distribution. Then, such xth quantile is denoted by Q GPD−NHPP (x). As mentioned previously, we assume that K 0 = Q GPD−NHPP (0.5), K 1 = Q GPD−NHPP (0.75) and K 2 = Q GPD−NHPP (0.95).
In this setting, based on Monte Carlo simulations, the cat bond price is estimated as 0.822337, so the relative percentage change (compared to the reference value from Example I.1) is equal to −0.9811134458 %. Then in this case, the change of the distribution of the single loss (but with properly estimated parameters) is less significant than, e.g., shift in the parameters discussed in Example I.3.

Example I.7
Similar to the previous cases, the relation among the cat bond price and the parameters ξ GPD , β GPD of the distribution of the single loss (given by GPD) is analyzed, so the prices for ξ GPD ∈ [0.5, 1.1] and β GPD ∈ [5.2 × 10 7 , 5.45 × 10 7 ] are found. The other parameters are the same as in Example I.6. Then once again, the parameters of the distribution of the  The relative percentage change of the cat bond price ( IB p (0)) and the relative percentage change of ξ GPD ( ξ GPD ) are compared to the reference cat bond price and ξ GPD from Example I.6 single loss are changed compared to the parameters fitted to the historical data. As it could be seen (see Fig. 5), the cuts of the graph seem to be hyperbolic (if β GPD is constant, the exact values are given in the third row of Table 6 for β GPD = 5.2 × 10 7 ) or almost linear (if ξ GPD is constant, see Table 7 for ξ GPD = 0.8) in the given intervals of values. The relative percentage changes of prices (compared to the reference value from Example I.6) for both parameters are considerable higher than in all the previous examples (more than 26 % relative reduction of the price if ξ GPD is decreased about 38 %). Because of the mentioned hyperbolic shape of function for the fixed β GPD , the possible estimation error or future changes of ξ GPD may lead to significant undervaluation of the cat bond price.

Cat bond pricing-fuzzy case
As noted in Sect. 6, also the fuzzy case of the cat bond prices may be interested for analysts and policyholders because of uncertainty caused by future behavior of the financial market. The relative percentage change of the cat bond price ( IB p (0)) and the relative percentage change of β GPD ( β GPD ) are compared to the reference cat bond price and β GPD from Example I.6 Therefore, the influence of some fuzzy parameters of the interest rate model on the price is more deeply analyzed using Algorithm 2 introduced in Sect. 7.2. The overall parameters of the model of losses are the same as in Example I.1-i.e., we assume NHPP of quantity of catastrophic events, and the value of each loss is modeled by lognormal distribution. The applied two-factor Vasicek model arises from parameters considered also in Sect. 8.1.1, and other assumptions about the cat bond are the same as in that setting.
The applied parameters are based on historical data, but because of the uncertainty related to future behavior, the experts' opinions should be incorporated into estimation of the prices [see, e.g., Nowak and Romaniuk (2014c)].
We focus on three important variables of the interest rate model, namely volatilities σ r , σ ε and the correlation coefficient ρ. They are important part of the model but they could be also wrongly estimated. In the following, only triangular fuzzy numbers are considered, but the similar numerical approach may be also used for other types of LR numbers.

Example II.1
For the parameters σ r , σ ε , ρ we use triangular fuzzy numbers related to the historical values from Puhle (2008). Then, we havẽ Monte Carlo simulations (see Sect. 7.2) lead to evaluation of the fuzzy cat bond price (see Fig. 6) which is almost triangular and symmetric in its shape. The analyst may be also interested in the behavior of fuzzy numbers. For example, instead of previously usedσ r , other symmetric triangular fuzzy values ofσ r may be used.  Fig. 8).
Then, the fuzzy cat bond price is almost symmetrical but clearly LR, not a triangular, number. Analysis similar to Example II.1 may be also conducted. The fuzzy cat bond prices evaluated for various asymmetrical values ofσ ε are illustrated in Fig. 9. Once again, for wider support ofσ ε , the output fuzzy price has also wider support.

Example II.3
For the fixedσ r = [0.0063, 0.0073, 0.0083] andσ ε = [0.0119, 0.0219, 0.0419], the cat bond prices for various types of triangular valuesρ are estimated. In the considered setting, the type of triangular numberρ (see Fig. 10) has rather moderate effect on the output and both symmetric or asymmetric values ofρ lead to similar cat bond prices.

Evaluation of the membership function
Now, we focus on numerical approaches considered in Sect. 7.3. In Sect. 8.3.1, bisection search is used to evaluate membership value α for the given current market price of the cat bondĈ t . Then, Algorithm 3 is applied in Sect. 8.3.2 to obtain membership functions for the extended advice choice function for the fixed actual cat bond price (i.e., if this current market price with the obtained fuzzy value is compared). Finally, Algorithm 4 is used in Sect. 8.3.3 to find the relevant membership functions for˜ for the whole interval of prices.  As explained in Sect. 7.3, special optimization procedure based on bisection search is used to find membership value α for the fixed current market price of the cat bondĈ t , i.e., to find μĨ B(t) Ĉ t solving the problem (OP2). We apply parameters from Sect. 8.2.1 to calculate α for variousĈ t . Some exemplary results for the fixed t = 0 are given in Table 8.

Example III.2
Using Algorithm 3, it is possible to evaluate the membership value α for the extended advice choice function˜ introduced in Sect. 6 for the single fixed current market cat bond price. In Table 9, there are values of α found for various current market pricesĈ 0 and possible decisions (to buy/to accumulate/to hold/to reduce/to sell) from the set V . Once again parameters from Sect. 8.2.1 are used to approximate the fuzzy cat bond price.

Example III.3
The whole shape of membership function for each decision from the set V could be also estimated and plotted using Algorithm 4. Applying the parameters from Sect. 8.2.1, the relevant membership functions for the decisions to buy/to accumulate and to hold/to reduce are illustrated in Figs. 11 and 12, respectively. Additionally, the financial analyst could take relevant decision based on possible recommendations for the fixed value α of the membership function. The example of such approach is presented in Fig. 13, where α = 0.95 is used. Therefore, knowing the current market priceĈ t , the analyst could take his decision based on similar sets.
If g(α) = Ĩ B (t) L α and h(α) = Ĩ B (t) U α are strictly monotonic functions (as in the previously considered examples), then the simplified approach could be used instead of Algorithm 4. In such case, Algorithm 2 gives the approximation of the fuzzy cat bond price which then leads to the output similar to the one summarized in Fig. 13

Conclusions
The increasing number of natural catastrophes leads to severe problems with stability of insurance industries. Therefore, new financial and insurance instruments are required, such as the catastrophe bonds.
In this paper, we extend our earlier results concerning cat bond pricing. As the model of spot interest rate, we apply the two-factor Vasicek model. We analyze the influence of various parameters on the cat bond price applying Monte Carlo methods.
Moreover, taking into account the uncertainty on the market, we derive fuzzy counterparts of the crisp pricing formulas and present an automated approach for decision making in fuzzy environment with illustrative examples.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.