Valuing catastrophe bonds involving correlation and CIR interest rate model

Natural catastrophes lead to problems of insurance and reinsurance industry. Classic insurance mechanisms are often inadequate for dealing with consequences of catastrophic events. Therefore, new ﬁnancial instruments, including catastrophe bonds (cat bonds), were developed. In this paper we price the catastrophe bonds with a generalized payoff structure, assuming that the bondholder’s payoff depends on an underlying asset driven by a stochastic jump-diffusion process. Simultaneously, the risk-free spot interest rate has also a stochastic form and is described by the multi-factor Cox–Ingersoll–Ross model. We assume the possibility of correlation between the Brownian part of the underlying asset and the components of the interest rate model. Using stochastic methods, we prove the valuation formula, which can be applied to the cat bonds with various payoff functions. We use adaptive Monte Carlo simulations to analyze the numerical properties of the obtained pricing formula for various settings, including some similar to the practical cases.


Introduction
Nowadays overwhelming risks caused by natural catastrophes, like hurricanes, floods and earthquakes, lead to severe problems of insurance and reinsurance industry.  Muermann 2008). Such extreme losses from a single catastrophic event cause problems relating to reserve adequacy or even lead to bankruptcy of insurers. For example, after Hurricane Andrew more than 60 insurance companies fell into insolvency (see Muermann 2008).
The main reason of the mentioned problems is related to assumptions used in classical insurance mechanisms. In traditional insurance models (see, e.g., Borch 1974) risk claims that are independent and small in relation to the value of the whole insurance portfolio (e.g. caused by car crashes) are the norm. Then the classic strategy of building portfolio (the grater the number of risks, the better quality of the whole portfolio) is justified by the law of large numbers and the central limit theorem (see, e.g., Borch 1974;Ermoliev et al. 2001). In the case of natural catastrophes, the sources of risks are strictly dependent on time and location. Additionally, problems with adverse selection, moral hazard and the cycles of prices of reinsurer's policies should be noted (see, e.g., Ermoliev et al. 2001;Finken and Laux 2009).
Therefore, new financial derivatives which connect both the financial markets and the insurance industry were developed. The main aim of these instruments is to securize the catastrophic losses, i.e. to transfer insurance risks into financial markets by "packaging" of risks into special tradable assets-catastrophic derivatives (see, e.g., Cummins et al. 2002;Freeman and Kunreuther 1997;Froot 2001;Harrington and Niehaus 2003;Nowak 1999;Nowak and Romaniuk 2010b, c, d;Nowak et al. 2012).
One of the most popular catastrophe-linked security is a catastrophe bond (known also as a cat bond or an Act-of-God bond (see, e.g., D'Arcy and France 1992;Ermolieva et al. 2007;George 1999;Nowak and Romaniuk 2009a;O'Brien 1997;Romaniuk and Ermolieva 2005;Vaugirard 2003). In 1993, catastrophe derivatives were introduced by the Chicago Board of Trade (CBoT). These financial derivatives were based on underlying indexes reflecting the insured property losses due to natural catastrophes reported by insurance and reinsurance companies. Then new approaches to development of the cat bonds were applied (see, e.g., Kwok 2008; Lee and Yu 2007).
The payoff received by the cat bondholder is linked to an additional random variable, which is called triggering point. This event (indemnity trigger, parametric trigger or index trigger) is usually related to occurrence of specified catastrophe (like hurricane) in given region and fixed time interval or it is connected with the value of issuer's actual losses from catastrophic event (like flood), losses modeled by special software based on the real parameters of a catastrophe, or the whole insurance industry index, or the real parameters of a catastrophe (e.g., earthquake magnitude or wind speeds in case of windstorms), or the hybrid index related to modeled losses (see, e.g., George 1999;Niedzielski 1997;Vaugirard 2003;Walker 1997). In the case of some cat bonds, the triggering point is related to the second or even the third event during a fixed period of time. Additionally, the structure of payments for the cat bonds depends also on some primary underlying asset (e.g. the LIBOR).
As noted by many authors (see, e.g., Ermoliev et al. 2001;Finken and Laux 2009;Vaugirard 2003), the cat bonds are important tools for insurers and reinsurers. Among other advantages, they stressed that using the cat bonds lowers the costs of reinsurance and reduces the risks caused by moral hazard.
The cash flows related to the cat bond are usually managed by special tailor-made fund, called a special-purpose vehicle (SPV) or a special purpose company (SPC) (see, e.g., Lee and Yu 2007;Vaugirard 2003). The hedger (e.g. insurer or reinsurer) pays an insurance premium in exchange for coverage in the case if triggering point occurs. The investors purchase an insurance-linked security for cash. The mentioned premium and cash flows are directed to SPV, which issues the catastrophe bonds. Usually, SPV purchases safe securities in order to satisfy future possible demands. Investors hold the issued assets whose coupons and/or principal depend on occurrence of the triggering point, e.g. the catastrophic event. If this event occurs during the specified period, the SPV compensates the insurer and the cash flows for investors are changed. Usually, these flows are lowered, i.e. there is full or partial forgiveness of the repayment of principal and/or interest. However, if the triggering point does not occur, the investors usually receive the full payment.
In the literature concerning the catastrophe bonds and their pricing many authors apply stochastic models. Among them one should mention two advanced approaches, where stochastic processes with discrete time are used:  within the framework of representative agent equilibrium and Reshetar (2008), where the payoff functions depend on catastrophic property losses and catastrophic mortality.
More authors apply stochastic models with continuous time. To incorporate various characteristics of the catastrophe process compound Poisson models are used in Baryshnikov et al. (1998). In this approach, no analytical pricing formula is obtained and the problem of change of probability measure in the arbitrage method is not discussed. However, advanced numerical simulations are conducted and analyzed. The authors of (Burnecki et al. 2003) correct the method proposed in Baryshnikov et al. (1998). In turn, the approach from Burnecki et al. (2003) is applied in Härdle and Lopez (2010) for the cat bonds connected with earthquakes in Mexico. In Albrecher et al. (2004) the doubly stochastic compound Poisson process is used and reporting lags of the occurred claims are incorporated to the model. The model behavior is analyzed with application of QMC algorithms. The arbitrage method for cat bonds pricing is used by Vaugirard (2003). He addresses the problem of non completeness of the market, caused by catastrophic risk, and non-traded insurance-linked underlyings in the Merton's manner (see Merton 1976). In the approach proposed in Lin et al. (2008) the Markov-modulated Poisson process is applied for description of the arrival rate of natural catastrophes. Jarrow in Jarrow (2010) obtained an analytically closed cat bond valuation formula, considering the LIBOR term structure of interest rates.
In Nowak and Romaniuk (2013a) we applied the approach similar to the one proposed in the Vaugirard's paper. We proved a generalized catastrophe bond pricing formula, assuming the one-factor stochastic diffusion form of the risk-free interest rate process. In contradistinction to the Vaugirard's approach, where catastrophe bonds payoffs were dependent on risk indexes, we considered the cat bond payoffs dependent only on the cumulated catastrophic losses. Moreover, we conducted Monte Carlo simulations to analyze the behavior of the valuation formula. The mentioned paper summarized and generalized our earlier results from Nowak and Romaniuk (2010b, c, d), Nowak et al. (2012). Shortly after our publication, results similar to ours were obtained in Ma and Ma (2013), where the authors assumed the one-factor Cox-Ingersoll-Ross (CIR) model of the risk-free interest rate.
In Romaniuk (2009b, 2013b) we considered the problem of cat bond pricing in fuzzy framework, incorporating uncertain financial market parameters to the model. Similar approach was also applied by us in Nowak and Romaniuk (2010a, 2013c, 2014, where the stochastic analysis, including the Jacod-Grigelionis characteristics (see, e.g., Shiryaev 1999;Nowak 2002), and the fuzzy sets theory were employed to find the European option pricing formulas.
In this paper we continue our considerations concerning valuation of the catastrophe bonds. We assume no arbitrage on the market and the possibility of replication of interest rate changes by financial instruments existing on the market. We use the martingale method of pricing. We apply the d-dimensional Brownian motion (with d ≥ 1) for description of the risk-free spot interest rate and the one-dimensional Brownian motion and the compound Poisson process to model an underlying asset I , connected with the cumulative catastrophic losses. We assume that I is similar to a synthetic insurance industry asset. Our contribution is threefold. First, we consider the catastrophe bonds with a generalized (in comparison to Vaugirard 2003; Nowak and Romaniuk 2013a) payoff structure, depending on I . Second, in contradistinction to our previous approaches, where the one-factor spot interest rate models were applied, the interest rate behavior is described by the multi-factor CIR model. Third, we assume the possibility of correlation between the Brownian part of the process I , used in the catastrophe bond payoff, and the components of the model of the interest rate. To our best knowledge such the approach assuming correlation structure has not been considered in the pricing literature.
There are several essential differences between the approach presented in this paper and the model of Vaugirard (2003). We apply the multi-factor Cox-Ingersol-Ross model of the risk-free spot interest rate, whereas in Vaugirard (2003) the one-factor Vasicek interest rate model is used. In Vaugirard's approach the catastrophe bond payoff depends on a physical risk index driven by the Poisson jump-diffusion process. Its Brownian part models the unanticipated instantaneous index change, reflecting causes that have marginal impact on the gauge. In turn, jumps, connecting with catastrophic events, increase the value of the risk index. In our approach we also use the jump-diffusion process to model an underlying asset, similar to a synthetic insurance industry asset. However, in contradistinction to Vaugirard (2003), its Brownian part plays a more important role, modelling, similarly as in Merton (1976), vibrations in price caused by temporary imbalance between supply and demand on the market. Jumps, connected with occurrences of catastrophic losses, decrease the value of an underlying instrument. For technical reasons we use a transformation of the process of an underlying asset for description of the bondholder payoff. The payoff structure in our model is much more general then the one considered in Vaugirard (2003). In particular, it is possible to use a wide class of functions for description of dependence between the bondholder payoff and the transformed underlying asset process. Finally, as we have mentioned above, in contradistinction to the Vaugirard's model and models of other authors, our approach enables taking into account the possibility of correlation between the Brownian part of the underlying asset and the Brownian motions modelling the interest rate behavior.
Since we use stochastic models of the spot interest rate and the underlying asset, stochastic analysis methods play the key role in derivation and proof of the cat bond pricing formula. In particular, the Girsanov theorem and Lévy's characterization of the Brownian motion is used. Furthermore, the correctness of the applied method of change of probability measure is proved in detailed way. The proposed by us payoff structure enables to use a wide class of functions describing dependence between the values of bondholder's payoff and the asset I , including stepwise, piecewise linear and piecewise quadratic one.
Apart from theoretical considerations, we conduct simulations to compare behavior of these models for different payoff structures. In numerical experiments we find the prices of the catastrophe bonds applying linear and quadratic payoff functions. To analyze the behavior of the obtained prices, we alter some parameters of the appropriate interest rate model and the model of value of catastrophic losses. This paper is organized as follows: Sect. 2 contains necessary notations and definitions concerning stochastic notions and processes used in the paper. Moreover, assumptions concerning the financial market are formulated. Section 3 contains generalized definition of the catastrophe bond payoff structure as well as description of the multi-factor CIR interest rate model. In Sect. 4 the catastrophe bond pricing formula is introduced and proved. Since the mentioned above risk-free interest rate is modeled by the multi-factor affine process, the underlying asset is defined by the stochastic jump-diffusion and the Brownian parts of both processes can be correlated, the derivation and proof of the valuation formula required application of stochastic analysis. Apart from Theorem 3, which is the main theorem of the paper, Lemma 2 is formulated to describe the catastrophe bond price at the moment zero and simplify computations of the expected bondholder's payoff. In Sect. 5 adaptive Monte Carlo simulations are conducted for the introduced formula of cat bond pricing. First, some necessary numerical algorithms are considered. Then the cat bond prices are estimated and analyzed for various settings, including the set of parameters similar to the practical case. Special attention is paid to the influence of the parameters of the underlying asset like correlation coefficients (which are important properties of the model considered in this paper) on the numerically evaluated price.

Stochastic and financial preliminaries
In this section we introduce some necessary notations, definitions and assumptions concerning stochastic models of catastrophe losses and financial market.
We denote by . the Euclidean norm in R d , i.e. for each vector x ∈ R d of the form Here denotes transposition so that x is a column vector. Moreover, we will use the symbol |||.||| to denote the Euclidean norm in R d+1 . R d×d denotes the space of d × d matrices of real numbers.
In the further part of the paper we will use the notion of quadratic covariance, which is generally defined for semimartingales (for details we refer the reader to Shiryaev 1999). Let T ⊆ [0, ∞) be a time interval.

Definition 1
We call a stochastic process X = (X t ) t∈T a semimartingale if it is representable as a sum where A is a process of bounded variation (over each finite interval [0, t]), M is a local martingale, both defined on a filtered probability space, satisfying the usual conditions. Definition 2 For two semimartingales X and Y , on a filtered probability space, the quadratic covariance process is the process where t 0 X s− dY s and t 0 Y s− dX s are stochastic integrals with respect to Y and X , respectively.
For description of losses caused by natural catastrophes and behavior of the risk-free spot interest rate on the market we apply stochastic processes with continuous time. In the paper we consider three different probability measures. For a probability measure M we denote by E M the expected value with respect to this measure. In particular, all the stochastic processes and random variables introduced in this section are defined with respect to a probability P.
All the economic activity will be assumed to take place on a finite horizon 0, T * , where T * is a positive constant. Let for a positive integer d be the standard d-dimensional Brownian motion. That is, each W i t is the one-dimensional Brownian motion and the different components W 1 t , W 2 t , . . . , W d t are independent. The process W X will be used for description of the interest rate on the market. We additionally consider the one-dimensional Brownian motion W I t t∈[0,T * ] , used in the further part of the paper for description of an underlying asset. We assume that for i = 1, 2, . . . , d W I and W i can be correlated with a correlation coefficient ρ i , i.e. the quadratic covariations have the form of independent and identically distributed non-negative random variables with finite expectation to describe values of losses during catastrophic events.
For each t ∈ 0, T * cumulative catastrophe losses until the moment t are modeled by the compound Poisson processÑ where N t is the standard Poisson process with a constant intensity κ > 0. Moments of jumps of the process N t correspond to moments of catastrophic events. We denote byN t the jump processN t =Ñ t − κe 1 t, where e 1 = 1 − E P e −U i . As we mentioned earlier, all the discussed above stochastic processes and random variables are defined on a probability space ( , F , P). We introduce the filtration are independent. Furthermore, the probability space with filtration , F , (F t ) t∈[0,T * ] , P satisfies the usual assumptions: the σ -algebra F is P-complete, the filtration (F t ) t∈[0,T * ] is right continuous and each F t contains all the P-null sets from F . By the symbol (B t ) t∈[0,T * ] we denote banking account satisfying the standard stochastic equation: is the risk-free spot interest rate, i.e. short-term rate for risk-free borrowing or lending at time t over the infinitesimal time interval [t, t + dt]. In the paper we assume that r is modeled by a time-homogenous d-dimensional Markov diffusion process X = (X 1 , X 2 , . . . , X d ) , given by the equation with the value space S ⊆ R d . The functions α : S → R d and σ : S → R d×d are sufficiently regular so that the above equation has a unique solution. We consider a particular diffusion model, i.e. the multi-factor Cox-Ingersoll-Ross model, where the spot risk-free interest rate process (r t ) t∈[0,T * ] has the form for an affine function r : S → R and X, defined in detail in Sect. 3.2. For all t ∈ 0, T * the banking account process B t has the form We assume that zero-coupon bonds are traded on the market and by the symbol B (t, T ) we denote the price at the time t, t ∈ 0, T * , of the zero-coupon bond with the face value equal to 1 and the maturity date T ≤ T * . Similarly as in Vaugirard (2003), we take into account the possibility of catastrophic events, using the jump-diffusion process Since our approach is general, we do not characterize precisely the instrument I . However, it can be interpreted as an instrument similar to a synthetic insurance industry underlying asset.
Moreover, we introduce the stochastic process which will be used in definition of the catastrophe bond payoff function in Sect. 3.1. We make the following assumptions concerning financial market: there is no possibility of arbitrage; there are no restrictions for borrowing and short selling; trading on the market takes place continuously in time; there are no transaction costs; lending and borrowing rates are equal and changes in the interest rate r can be replicated by existing financial instruments.

Payoff structure
The payoff structure of the catastrophe bond is described by classes W, and K defined below.
We fix a positive integer n ≥ 1, a face value of the catastrophe bond Fv > 0 and a maturity date of the cat bond T ∈ [0, T * ].
The class of sequences where 0 ≤ w 1 , w 2 , . . . , w n and n i=1 w i ≤ 1, is denoted by W. The partial sums of w ∈ W are denoted by is the class of sequences of functions ϕ = (ϕ 1 , ϕ 2 , . . . , ϕ n ) fulfilling the following conditions: In particular, we consider the following subclasses of : Clearly, 1,l and 1,q are subclasses of 1 .
Finally, K is the class of increasing sequences We proceed to define the catastrophe bond payoff function. Let w ∈ W, ϕ ∈ and K ∈ W. We introduce an auxiliary function satisfying the following assumptions: Definition 3 Let w ∈ W, ϕ ∈ and K ∈ W. We denote by I B (w, ϕ, K ) the catastrophe bond with the face value Fv, the maturity and the payoff date T if its payoff function is the random variable ν w,ϕ,K given by the equality The payoff function ν w,ϕ,K of I B (w, ϕ, K ) will be called stepwise (piecewise linear or piecewise quadratic) if ϕ ∈ 0 (ϕ ∈ 1,l or ϕ ∈ 1,q ).
The following remark shows basic facts concerning the cat bond defined above. The presented formulas are obtained by straightforward computations.
Remark 1 The catastrophe bond I B (w, ϕ, K ) has the following properties: 1. The general formula describing the payoff as a function ofĪ T can be written in the form In particular, 2. IfĪ T is relatively small (i.e.,Ī T ≤ K 0 ), the bondholder receives the payoff equal to its face value Fv. 3. IfĪ T > K n , the bondholder receives the payoff equal to Fv(1 − w (n) ).
. . , n, the bondholder receives the payoff equal to In case of the stepwise payoff function this payoff is constant and equal to Fv

The multi-factor Cox-Ingersoll-Ross interest rate model
Multi-factor affine interest rate models were introduced by Duffie and Kan (1996). Their paper (see Duffie and Kan 1996) is regarded as a cornerstone in the interest rates term structure theory. Dai and Singelton (2000) provided classification of the multi-factor affine interest rate models and reasoning on their structure. The popularity of the mentioned models follows from their tractability for bond prices and bond option prices. A multi-factor affine model of the interest rate is described by a time homogeneous diffusion model given by where κ and are constant d × d matrices, ϕ = (ϕ 1 , ϕ 2 , . . . , ϕ d ) is a constant vector, for i = 1, 2, . . . , d υ i are constants and are constant vectors, We also assume that there is a real constant ξ 0 and a constant vector As we mentioned earlier, in this paper we consider the multi-factor Cox-Ingersoll-Ross model, which is an affine interest rate model of the form (6) with ϕ i > 0, ξ i = 1 and υ i = 0 for i = 1, 2, . . . , d as well as

Catastrophe bond pricing formula
Our aim in this section was to prove the catastrophe bond pricing formula. We will apply the following version of the Levy theorems (see, e.g., Ikeda and Watanabe 1989;Shreve 2004).

Theorem 1 Let M be a martingale relative to a filtration
Then M is the Brownian motion.
In Vaugirard (2003) the author considered a simple form of the catastrophe bond payoff function. The triggering point was defined as the first passage time through a level of losses K of a natural risk index I . He assumed that if the triggering point does not occur, the bondholder is paid the face value Fv; and if the triggering point occurs, the payoff is equal to the face value minus a coefficient in percentage w, i.e. Fv(1 − w). Bondholders were regarded to be in a short position on a one-touch up-and-in digital option on I and, similarly as in case of options, the martingale method was used to find the catastrophe bonds valuation expression. In our approach we also use the conditional expectation with respect to equivalent risk-neutral measure to obtain the analytical form of the cat bond pricing formula. According to our earlier definitions, the model considered by us has the more general payoff structure, the underlying asset I is connected with the insurance industry and there is the possibility of correlation between the continuous parts of the processes describing r and I . Now we formulate and prove the main theorem concerning catastrophe bond pricing.
Theorem 3 Let (r t ) t∈[0,T * ] be a risk-free spot interest rate process given by the multi-factor CIR model r t = r (X t ), t ∈ 0, T * , where X is the vector process with parameters described in the previous section. Let w ∈ W, ϕ ∈ , K ∈ K, T ≤ T * and let I B w,ϕ,K (0) be the price at time 0 of I B (w, ϕ, K ). Then there exists a probability measure Q F , equivalent to P, such that where where the vector process X is described by the stochastic equation with the matrixκ of the form In formulas (11) and (12) Proof From the theory of assets pricing it follows that for risk-neutral equivalent probability measure Q. Let W d+1 be the Brownian motion obtained by Lemma 1 for k = d, Z i = W i , 1 ≤ i ≤ d, and Z d+1 = W I . Then W d+1 satisfies the equality: The change of probability measure is described by the Radon-Nikodym derivative dQ d P = Z T * P-a.s., where for the market price of risk processesλ X t andλ d+1 t of the form In formula (16)λ = λ 1 ,λ 2 , . . . ,λ d is a constant vector and then (Z t ) t∈[0,T * ] is a martingale with respect to P and Q is a probability measure equivalent to P. Let us assume that the equality (17) is satisfied. Then from the Girsanov theorem (see, e.g., Karatzas and Shreve 1988, Chapter 3.5.A.) it follows that there exist the two independent Q-Brownian motions: d-dimensionalW X = W 1 ,W 2 , . . . ,W d and one-dimensionalW d+1 , satisfying the equalities: SinceW I t is a continuous martingale starting from 0 and W I ,W I t = t, by Theorem 1, it is the Q-Brownian motion. Moreover, for t ∈ [0, T * ] and 1 ≤ i ≤ d, W I ,W i t = ρ i t. From (14) it follows that Under Q the vector process X t is given by the equation (19) implies and after reformulation It remains to prove (17). To this end, we apply an idea similar to the one used in Cheridito et al. (2007).
be the solution of the equation with respect to P. Let Let n ≥ 1. We introduce the following stopping times: The process Therefore, the probability measure Q n , defined by the Radon-Nikodym derivative dQ n d P = Z n T * , where is a probability measure equivalent to P and the process W n,X = W n,1 ,W n,2 , . . . ,W n,d , satisfying the equality dW n,X t = dW X t +λ n,X t dt, is the Q n -Brownian motion. Furthermore, the process with respect to Q n has the same distribution as the process X t∧τ n t∈[0,T * ] with respect to P. Therefore, Q n τ n = T * = P τ n = T * .
Moreover, one can check that Applying the monotone convergence theorem (see, e.g., Billingsley 1986, Theorem 16.2), we obtain the equality Since, for each n ≥ 1, E P Z T * I {τ n =T * } = Q n τ n = T * , the equalities (24), (25) imply (17). For t = 0 the zero-coupon bond price has the form Moreover, from Munk (2011) it follows that B (t, T ), t ∈ [0, T ], satisfies the equation where and its solution has the form (10), which finishes the proof of the assertion (i). We introduce the next probability measure Q F , equivalent to Q, given by the following Radon-Nikodym derivative: where are the independent Q F Brownian motions. Sincẽ is a continuous martingale starting from 0 and W I ,W I t = t for t ∈ [0, T * ], Theorem 1 implies thatW I t is the Q F -Brownian motion. Moreover, for i = 1, 2, . . . , d and t ∈ [0, T * ], W I ,W i t = ρ i t. Equalities (18), (26) and (27) imply the equality dW I t = dW I t − ρ σ T t dt. Therefore, the processes (20) and (21) take the following form with respect to Q F : This finishes the proof of the assertion (ii). Clearly, and application of formula (29) to (28) gives (9).
By straightforward computations, applying Theorem 3, we obtain the following lemma concerning a detailed form of the catastrophe bond price at the moment 0.
The above lemma simplifies the numerical computations of the catastrophe bond price.

Numerical simulations
In order to analyze the behavior of prices of the cat bonds, the Monte Carlo simulations are conducted in this section. To utilize the obtained general pricing formula (9) proved in Theorem 3, iterative schemes of simulations for the process I t given by (11) and the process X t given by (12) are applied with fixed time step t. Because the final estimator of the price especially depends on supremum of the generated trajectory of the processĪ t defined by (4), the adaptive approach with adjustment of the length of t was introduced. The necessary algorithms are considered in Sect. 5.1. In the following we illustrate the possibility of pricing cat bonds in various parametric settings via numerical computations, despite the complex nature of the formulas considered in Theorem 3 and Lemma 2. We start our considerations from the simplified, synthetic setup discussed in Sect. 5.2. In some parts this first setting is similar in nature to the one considered in Vaugirard (2003) or other financial papers (see, e.g., Nowak and Romaniuk 2010a). The following discussion enables us to track down the most important details of behavior of prices of the considered types of the catastrophe bonds.
Then in Sect. 5.3 we turn to other approach, which is more complex, real-life setup, because it is partially based on the model of interest rates and the model of catastrophic events adapted from Chen and Scott (2003) and Chernobai et al. (2006), respectively. In these two papers the real-life data were considered and the parameters of the related statistical models were estimated for this data. Therefore, we show that even for more complex setting which is closer to the practical cases, the evaluation and analysis of the prices of the considered cat bonds is possible and leads to important conclusions for practitioners. Also the structure of payments of the cat bond is statistically analyzed in this case.
In the following, one-or multifactor CIR models of interest rates are discussed. We also assume that the generated losses U i are of a catastrophic nature-they are rare, but each loss has a high value. Therefore, the quantity of losses is modeled by HPP (homogeneous Poisson process) and the value of each loss is given by a random variable with a relatively high expected value and variance (i.e., high risk with high variability). We limit our considerations to the case when the value of each loss is modeled by lognormal distribution. This distribution is commonly used in simulations of risk events in insurance industry. However, other distributions, e.g. Weibull, gammma, GEV (see Chernobai et al. 2005;Furman 2008;Hewitt and Lefkowitz 1979;Hogg and Klugman 1983;Melnick and Tenenbein 2000;Papush and Patrik 2001;Rioux and Klugman 2006), or simulations based on historical records (see Ermolieva and Ermoliev 2005;Pekárová et al. 2005) are possible and they can be easily incorporated into the approach presented in this paper.
We assume that the face value of the bond in each numerical experiment is set to 1 (i.e. one monetary unit assumption is used) and the trading horizon of the catastrophe bond is set to 1 year. The starting value I 0 of the process (3) is equal to 1. In each experiment we generate N = 1,000,000 simulations. The set of other necessary parameters of the catastrophe bond, the model of interest rates and the model of catastrophe events are described in details for each analysis.
In our considerations we focus on the piecewise linear payment function 1,l of the cat bond as defined in Remark 1. As previously noted, the price of the cat bond for this type of the payment function is directly related to supremum evaluated for the processĪ t . Therefore, the time step t applied in the iterative Monte Carlo scheme is adapted according to this value. Usually, t = t long = 0.02 is set which is close to 1-week cycle (for assumed T = 1). But if the value of the processĪ t is inside the interval of the values close to the first triggering point K 0 or the last one K n , the time step is shortened to t = t short = 0.005. Therefore, the obtained estimator has better quality and the whole numerical procedure is more flexible. Of course, additional moments related to the jumps caused by the catastrophic events U i are also taken into account in our approach as explained further.

Algorithms
As indicated by Theorem 3 and Lemma 2, the evaluated cat bond price I B (w, ϕ, K ) depends on three processes: the jump processÑ t defined by (2) (or its transformationN t , equivalently), the price of the underlying asset I t with respect to the probability measure Q F described by (11) and the interest rate model X t given by (12). The processes I t and X t are correlated via the Brownian motionsW X t andW I as defined by (13). Therefore, to apply the pricing formula (9), two types of simulations should be used. During the first one (Algorithm 1), the trajectory ofÑ t is generated. During the second one (Algorithm 2), both of the trajectories I t and X t are generated jointly using some input from the first type of simulation. Then, based on all of the trajectories, the expected value E Q F ν w,ϕ,K given by (30) is estimated via Monte Carlo approach (Algorithm 3). Additionally, the discounting factor B (0, T ) given by (10) is evaluated. Eventually, the main formula (9) could be used, merging these two outputs.
We start from description of the algorithm which is used to generate the processÑ t . Because HPP is applied, then the intervals between the consecutive jumps are given by iid random variables from exponential distribution with the parameter κ (see, e.g., Romaniuk and Nowak 2015). To generate the jumps U i , some fixed random distribution should be used. In the case of the lognormal distribution considered in this paper the relevant algorithms are widely studied in the literature (see, e.g., Romaniuk and Nowak 2015). Then we have the following steps which generate single trajectory of the jump process:
Step 2 Generate s from the exponential distribution with the parameter κ.
Step 4 Generate U from the distribution of losses. Set j = j + 1,Ñ t j =Ñ t j−1 + U, t j = t j−1 + s. Storẽ N t j , put t j into the sequence t jumps in increased order.
Step 5 Return to the step 2. Output The trajectory ofÑ t and the jump moments t jumps .
The transformation ofÑ t into the processN t is straightforward if the expected value E P e −U for the distribution of the jump U is at least numerically known.
In the model considered in this paper there is embedded dependency between the processes I t and X t . From (13), it is related to the correlation matrix of d + 1 dimensional normal distribution given by ⎛ Using Cholesky decomposition the relevant trajectories ofW 1 t ,W 2 t , . . . ,W d t andW I t could be simulated (see, e.g., Romaniuk and Nowak 2015) for the given set of times. Because of the special form of (33),W I t is generated as a linear combination of independent normal random variables used in simulation ofW 1 t ,W 2 t , . . . ,W d t and one additional normal sample.
In order to simulate the trajectories I t and X t , we use the iterative scheme for 0 = t 0 < t 1 < · · · . First of all, the time step t = t j+1 − t j depends on the processĪ t in which supremum of I t is used. If for some t the value of I t is such that I 0 I t is inside the interval [K 0 −ε, K n +ε] (for the fixed parameter ε > 0), then the time step is shortened to t short > 0; Otherwise, it is set to t long > t short . Such approach improves the efficiency of numerical simulations. Additionally, into the set of times for which I t and X t are generated, the sequence of moments of jumps t jumps from Algorithm 1 should be also incorporated.
In the considered setup Euler schemes are then used. From (12) and the assumptions about (7) introduced in Sect. 3.2 we get where i = 1, . . . , d and W 1 , . . . , W d are iid standard normal variables. From (11) and the assumptions introduced in Section 2 we have In the above formula, W I t j is increment of the Brownian motionW I t generated using the mentioned earlier linear combination of W 1 , . . . , W d with one additional normal variable W d+1 , e.g. if d = 2 then we have ρ 1 W 1 + ρ 2 W 2 + 1 − ρ 2 1 − ρ 2 2 W 3 . Moreover, where Ñ t j is increment of the jump processÑ t obtained using Algorithm 1. These considerations lead us to the following steps:

Algorithm 2
Input Parameters of: the multi-factor CIR model, the Poisson process (κ), the distribution of losses, the catastrophe bond (K 1 , . . . , K n ), the Brownian motions; the starting value I 0 , the accuracy rule (ε).
Step 1 Run Algorithm 1 to obtain the trajectory ofÑ t and its jump moments t jumps .
Step 3 If Step 4 Let t j+1 = t j + t and put t j+1 into the sequence of moments t all in increased order.
Step 5 Remove the earliest moment t j+1 from the sequence t all . If t j+1 > T , then t j+1 = T . Let t = t j+1 − t j .
Step 6 Find N t j using (36) and data from the step 1.
Step 8 Evaluate I t j using (35). Let I t j+1 = I t j I t j . Store I t j+1 .
Step 9 If t j+1 = T , then return the obtained trajectories I t and X t .
Step 10 Let t j = t j+1 , j = j + 1. Return to the step 3. Output The trajectory of I t and X t .
Then the sampled trajectory of I t is transformed intoĪ t . It allows us to obtain values necessary for evaluation of (30) (e.g. I K i−1 <Ī T ≤K i ). The last phase consists of approximation of the expected value E Q F ν w,ϕ,K , considered in Lemma 2 via crude Monte Carlo estimator (see, e.g., Romaniuk and Nowak 2015), and application of the main formula (9).

Algorithm 3
Input Parameters of: the multi-factor CIR model, the Poisson process (κ), the distribution of losses, the catastrophe bond, the Brownian motions; the starting value I 0 , the accuracy rule (ε), the number of simulations n.
Step 1 Run n times Algorithm 1 and Algorithm 2 to sample trajectoriesĪ Step 2 Find estimators of ψ i (see (31)) and e i (see (32)) using relevant averages based on the samplē I (1) Step 3 Evaluate B (0, T ) (given by (10)) and approximate E Q F ν w,ϕ,K (see (30)) using the estimators found in the step 2.
Step 4 Find the cat bond price with the pricing formula (9). Output The cat bond price I B w,ϕ,K (0).

Simplified setup
In Vaugirard (2003) the catastrophe derivatives written on the catastrophe index were considered. As the model of interest rates, Vasicek model was used. The intensity of catastrophic events was modeled by HPP, and the lognormal distribution described the value of single catastrophic loss. Then the simplified setup with intuitive parameters for the mentioned models was discussed to analyze the behavior of the cat bond prices.
In this paper, the process I t for the instrument similar to a synthetic insurance industry underlying asset with its transformationĪ t defined by (4) is considered. Therefore, we start our numerical analysis also from the simplified setup, which is close in its nature to Vaugirard (2003), to emphasize the most important features in the evaluation of the cat bond prices.
Model I: The relevant parameters of this pricing model are enumerated in Table 1. Similarly to Vaugirard (2003), in this simplified setup the intuitive values for the interest rate  (2003). The parameters of the lognormal distribution and the intensity of HPP are the same as in some of the analysis in Vaugirard (2003). Also the simple payment function with only two triggering points is considered. The value of reduction coefficient of the payoff w 1 = 0.9 is very high to emphasize the reduction of payment of the cat bond if the triggering point occurs. Such value is also similar to the one used in Vaugirard (2003). For illustrative purposes, two straightforward values of the triggering points K 0 and K 1 are also set (see Table 1).
The parameters of the Brownian component of the process I t are also presented in Table 1. To take into account the dependency between the behavior of the trajectory of the underlying instrument and the interest rates as described by Theorem 3, the correlation coefficient ρ 1 between the Brownian motions of these processes is set to 0.5. The variability σ I = 0.2 is relatively high comparing to the volatility of the interest rates and the parameters of the process of the catastrophic events.
Then the estimated price of the catastrophe bond in this case is equal to 0.767317. Model I, Analysis I The intensity of HPP is important parameter of the model of catastrophic events. Therefore, in Vaugirard (2003) the dependency between the price of the catastrophe bond and the intensity κ HPP is analyzed for a few simple cases. We also adopt the similar approach but conduct the relevant numerical simulations for the whole interval κ HPP ∈ [0.4, 1.6] instead of a few values as in Vaugirard (2003) (see Fig. 1 for the graph of the obtained cat bond prices). The other parameters in this analysis are the same as in Table  1. As it could be seen, the cat bond price is strictly decreasing function of κ HPP .
Model I, Analysis II In the model of the process I t , especially two parameters are important comparing to approaches considered in other papers: the volatility of the Brownian motion of the underlying asset σ I and the correlation coefficient ρ 1 between the Brownian motion W I t and the behavior of the model of the interest rate. Therefore, the influence of these parameters on the obtained estimator of the cat bond price should be analyzed.
The cat bond prices obtained for the wide range of the mentioned parameters, i.e. σ I ∈ [0.1, 0.9] and ρ 1 ∈ [0.1, 0.9] may be found in Fig. 2. The other parameters are the same as in Table 1. As easily seen, the price is the decreasing function of σ I . However, the influence of ρ 1 is not so straightforwardly noticeable from Fig. 2. Only using the single cut of the relevant surface from Fig. 2 for the fixed value σ I = 0.6 with new, appropriate scale the dependency between the parameter ρ 1 and the cat bond price is easier to found (see Fig. 3). Then the price is also decreasing function of ρ 1 , but the influence of σ I on the obtained estimator is more significant in this setting.
Model I, Analysis III Usually, the estimation of the distribution of the single catastrophic event is based on historical data. It is possible that the future jumps in the relevant process will follow other patterns or there could be some error in estimation procedure. Therefore,  Fig. 4). The other parameters are the same as in Table  1. As it may be seen, the cat bond price is decreasing function of both μ LN and σ LN . Model I, Analysis IV The enterprise which issues the catastrophe bond may be also interested in the dependency between the price of such instrument and the parameters of its payment function. For example, behavior of the cat bond price for various values of the coefficient w 1 may be analyzed. Example of the output of the related simulations can be  Table 1) the cat bond price is almost linearly decreasing function of w 1 . Model I, Analysis V As it may be seen from the formula (11), there is important relation between the processes I t and X t . Therefore, the parameters of the model of interest rate also affects the final price of the catastrophe bond. For example, such influence may be seen if the price is numerically evaluated for various values of (see Fig. 6). The other parameters are the same as in Table 1. Then for the given interval ∈ [0.1, 1.0] the estimated price is explicitly non-linear convex function.

Complex setup
The cat bond pricing is also possible for the more complex setup which is closer to the reallife cases. Therefore, the two-factor CIR model of interest rates is applied. The parameters of this model were estimated in Chen and Scott (2003) based on monthly data of the Treasury bond market using Kalman filter. Also the parameters of the Poisson process and the applied distribution of the value of the single loss are based on real-life data in this setting. These parameters are adapted from Chernobai et al. (2006), where the information of catastrophe losses in the United States provided by the Property Claim Services (PCS) of the ISO (Insurance Service Office Inc.) and the relevant estimation procedure for this data are considered. Triggering points K 0 = Q(0.5), K 1 = Q(0.75), Values of losses coefficients For each catastrophe, the PCS loss estimate represents anticipated industrywide insurance payments for different property lines of insurance covering. An event is noted as a catastrophe when claims are expected to reach a certain dollar threshold. Initially, this threshold was equal to $ 5 million; then due to economic reasons it was increased. The PCS loss index is close in its nature to the "catastrophic part" given byÑ t for the process of the underlying asset considered in our paper. Additionally, the PCS loss process has important practical significance, because it is used as triggering point in many financial derivatives (see, e.g., Monti and Tagliapietra 2009;Schradin 1996).
Model II As previously noted, the two-factor CIR model based on the parameters from Chen and Scott (2003) is considered (see Table 2). The value of each catastrophic loss U i is modeled by the lognormal distribution and the number of losses N t is described by HPP. These parameters are adapted from the values estimated in Chernobai et al. (2006).
In this setting the triggering points for the payment function are connected with exceeding the limits given by quantiles Q(x) of the mentioned "catastrophic part" of the process I t , namely exp Ñ T , which could be derived from the formula (4) for t = T . Then Q(x) is x-th quantile for the random variable exp Ñ T if the value of each loss is described by the lognormal distribution and the number of events is given by the homogeneous Poisson process. Four different triggering point with three values of payment coefficients w i are used (see Table 2). The behavior of such complicated payment function is also possible to analyze using numerical simulations which is done in Model II, Analysis I.
The Brownian motions of the two-factor CIR interest rate model depend on the related part of the process I t via the correlation coefficients ρ 1 and ρ 2 . Because the values of these correlation coefficients are statistically low (see Table 2); therefore, the relevant influence is not significant in this setting.
Then the evaluated cat bond price is equal to 0.770507. Model II, Analysis I Apart from the price of the cat bond, the structure of its possible payments is interesting for the issuer of such financial instrument. Using Monte Carlo simulations these payments can be analyzed, e.g. histogram (see Fig. 7), box-and-whisker plot (see Fig. 8) and descriptive statistics (see Table 3) for the cat bond price described by the parameters from Table 2 can be directly obtained. As it may be seen from the mentioned histogram, almost 40 % of frequency of final payments is equal to the face value of the cat  bond. This value also constitutes mode. On the other hand, about 6 % of payments is close to zero, i.e. the trajectory of the processĪ t exceeds or is very close to the last triggering point. The rest of the payments are very linear in their behavior except for the values related to 1 − w 1 = 0.7 and 1 − w 1 − w 2 = 0.4. The histogram and the box-and-whisker plot are highly left-skewed. Model II, Analysis II As previously noted, the parameters of the Brownian motion used in the process I t are significant part of the model considered in this paper. Therefore, also for the more complex and more "real-life" setting described by Model II, the dependency among these parameters and the cat bond price should be analyzed.
As it may be seen from Fig. 9, the cat bond price is decreasing function of the correlation coefficient ρ 1 , if ρ 2 = 0.2 and σ I = 2 are set. This behavior is similar to the one noted for Model I, Analysis II. However, the dependency between the cat bond price and the second correlation coefficient ρ 2 , if ρ 1 = 0.2, is not so straightforward (see Fig. 10). Because of the symmetrical nature of the two-factor CIR model, this situation seems to be related to lower values of the parameters of the second term of the model of interest rates, i.e. ϕ 2 , κ 2 , 2 , considered in this setting. Model II, Analysis III Apart from correlation coefficients, the volatility σ I is important part of the considered process I t . Therefore, as in Model I, Analysis II, the dependency between the cat bond prices and this parameter is considered. The relevant estimators of the prices of cat bond for σ I ∈ [0.1, 1.1] are plotted in Fig. 11. In this setting the correlation coefficients have slightly higher values than in previous analysis and are set to ρ 1 = 0.4 and ρ 2 = 0.4. The other parameters are the same as in Table 2. Then the cat bond prices are decreasing concave function of σ I .

Conclusions
Increasing number of natural catastrophes leads to problems of insurance and reinsurance industry. Since classic insurance mechanisms could be inadequate for dealing with consequences of extreme catastrophic events, even a single catastrophe could result in bankruptcy or insolvency of insurers and reinsurers. Therefore, new financial instruments were introduced to transfer risks from insurance to financial market. The catastrophe bond is an example of such financial instruments. This paper is devoted to pricing the catastrophe bonds with a generalized payoff structure. We assume that the bondholder's payoff depends on an underlying asset described by a Levy jump-diffusion. The risk-free spot interest rate is driven by the multi-factor Cox-Ingersoll-Ross model. We take into account the possibility of correlation between the Brownian part of the log-price of the underlying asset and the components of the CIR interest rate model. Applying methods of stochastic analysis, we proved the general catastrophe bond valuation expression, which is formulated in Theorem 3. The obtained pricing formula can be used for the cat bonds with various payoff functions. In particular, stepwise, piecewise linear or piecewise quadratic payoff functions are special cases of the considered payoff structure. We also formulated Lemma 2 to describe in a more detailed way the catastrophe bond price at the moment zero and simplify computations of the expected bondholder's payoff.
The considered pricing formulas are then used in the adaptive, iterative Monte Carlo simulations to analyze some numerical properties. The prices are estimated for various settings-the more simplified setup, which is similar to the one considered in Vaugirard (2003), and the more complex case, which is closer to the practical applications and is based on real-life parameters adapted from Chen and Scott (2003) and Chernobai et al. (2006). In our analysis of the prices various parameters are taken into account, like parameters of the single loss distribution, the intensity of the number of the catastrophic events, the shape of the payment function, etc. During estimation of the cat bond prices special attention is paid to variables which are specific features of the model of the process considered in this paper like the correlation coefficients and the volatility of the Brownian motion. The presented analysis may be helpful for experts involved in construction phase of the new catastrophe bond or when the financial instrument is already available on the market.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.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.