# Effect of earthquake sequences on risk-based catastrophe bond pricing

Harsh Mistry, Andres Hernandez, Philippe Guéguen, Domenico Lombardi

► To cite this version:

Harsh Mistry, Andres Hernandez, Philippe Guéguen, Domenico Lombardi. Effect of earthquake sequences on risk-based catastrophe bond pricing. Risk Analysis, 2024, 10.1111/risa.14288. hal-04500791

HAL Id: hal-04500791

https://hal.science/hal-04500791v1

Submitted on 14 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Original Article

# Effect of earthquake sequences on risk-based catastrophe bond pricing

Harsh K. Mistry<fnref n="1" />

<figure><img src="image_2.png" /></figure>

Andres Hernandez<fnref n="2" />

Philippe Guéguen<fnref n="2" />

Domenico Lombardi<fnref n="1" />,<fnref n="3" />

Correspondence

Harsh K. Mistry, Department of Mechanical, Aerospace and Civil Engineering, University of Manchester, Manchester, UK.
Email: harsh.mistry@manchester.ac.uk, hmistry@verisk.com, harshmistry96@gmail.com

Funding information

University of Manchester

## Abstract

Catastrophe bonds (cat bond in short) are an alternative risk-transfer instrument used to transfer peril-specific financial risk from governments, financial institutions, or (re)insurers, to the capital market. Current approaches for cat bond pricing are calibrated on seismic mainshocks, and thus do not account for potential effects induced by earthquake sequences. This simplifying assumption implies that damage arises from mainshocks only, while aftershocks yield no damage. Postearthquake field surveys reveal that this assumption is inaccurate. For example, in the 2011 Christchurch Earthquake sequence and 2016–2017 Central Italy Earthquake sequence, aftershocks were responsible for higher economic losses when compared to those caused by mainshocks. This article proposes a time-dependent aggregate loss model that takes into account seismicity clustering and damage accumulation effects in the computation of damage. The model is calibrated on the seismic events recorded during the recent 2016–2017 Central Italy Earthquake sequence. Furthermore, the effects of earthquake sequence on cat bond pricing is explored by implementing the proposed model on five Italian municipalities. The investigation showed that neglecting time-dependency may lead to higher difference (up to 45%) in the cat bond price when compared to standard approaches.

## KEYWORDS

average annual loss, catastrophe bonds, earthquake sequences, Markov chain Monte Carlo, time-dependency, seismic risk assessment

# 1 | INTRODUCTION

Earthquakes have impacted communities around the globe, resulting in significant death tolls and socioeconomic losses. With 18,426 casualties and over US$250 billion of economic losses, the 2011 Tohoku Earthquake is the deadliest earthquake on record (Swiss Re, 2023). Such impact on loss of life can be minimized by adopting earthquake-resistant designs and early warning systems (Allen & Melgar, 2019; Cochran et al., 2018). On the other hand, economic losses can be managed by using financial risk-transfer mechanisms such as earthquake insurance and catastrophe bonds (cat bond in short). Earthquake insurance is mostly used by private individuals and asset owners to protect their assets against financial damage caused by earthquakes in exchange of a regular and known payment called premium. On the other hand,

cat bonds are used by insurance companies, governments, or financial institutions to transfer part, or the full financial risk, to the capital market. The investors purchase cat bonds thanks to their relatively high rates of return and low time to maturity (typically, 3 to 5 years). Figure 1 shows the economic impact caused by some of the recent earthquakes between 2009 and 2020. Except for the 2011 Christchurch Earthquake, the insurance penetration is below 33%, resulting in that higher proportion of uninsured losses for most seismic events. Some of the factors that govern such low insurance penetration are uncertainties when pricing insurance policies and cat bonds, low risk perception, and household and income demographics and affordability (Kajwang, 2020; Kelly et al., 2020; Pothon et al., 2019). Among these factors, pricing of cat bonds depends on the modeled loss data. Studies related to pricing cat bonds mainly focused on either using global

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2024 The Authors. *Risk Analysis* published by Wiley Periodicals LLC on behalf of Society for Risk Analysis.

<footnote n="1">¹Department of Mechanical, Aerospace and Civil Engineering, University of Manchester, Manchester, UK</footnote>

<footnote n="2">²ISTerre, Université Grenoble Alpes, Université Savoie Mont-Blanc, CNRS, IRD, Université Gustave Eiffel, Grenoble, France</footnote>

<footnote n="3">³Disaster Resilience Solutions Ltd., London, UK</footnote>
<figure><img src="image_3.png" /><figcaption>FIGURE 1 Economic loss for selected earthquakes around the world between 2009 and 2022.</figcaption></figure>

historical loss data or simulated loss data computed from catastrophe models.

Generally, cat bonds are priced using a stochastic process assuming arbitrage-free market and continuous time conditions. Such process follows a method used in pricing credit derivatives. Some of the early studies on cat bonds have focused on pricing zero-coupon and coupon cat bonds using 10-year global historical loss data (Burnecki & Kukla, 2003; Ma & Ma, 2013). Cardenas et al. (2007) calibrated cat bond for Mexican earthquake using parametric type of trigger mechanisms, which rely on the characteristics of the earthquake in a particular zone (also called cat-in-a-box). The main drawback of this approach is that it introduces higher levels of errors. Later these errors were minimized by using optimization algorithms that account for correlation between trigger parameters and event losses (Bayliss et al., 2020; Franco, 2010). Härdle and Cabrera (2010) introduced a hybrid type of trigger mechanism and calibrated earthquake cat bond price for Mexico. This approach is based on the combination of modeled earthquake losses and parametric triggers based on parameters such as moment magnitude exceeding a certain threshold. Shao et al. (2016), Shao et al. (2017) proposed a framework for pricing nuclear catastrophe risk bonds with three layers of risk categories: major, accident and incident; the transition from one damage state to another was modeled using a semi-Markov chain Monte Carlo approach. Hofer et al. (2019) introduced risk-based cat bond pricing formulations to account for uncertainties in the event arrival rate and loss severity. Based on these formulations, Hofer et al. (2020) presented a study where cat bonds were calibrated using simulated loss data obtained from simplified cat models. Later, Mistry and Lombardi (2022) proposed further enhancement in the risk-based cat bond pricing pro-

cedure by developing high-resolution hazard and exposure models. Recently, Mistry and Lombardi (2023) proposed a stochastic approach to account for uncertainty associated to asset location and attributes within the risk-based cat bond pricing procedures.

So far studies on risk-based cat bond pricing have used time-independent catastrophe models, which are based on the assumption that the occurrence of events within a given year is independent of each other (Hofer et al., 2020; Mistry & Lombardi, 2022, 2023). Such assumption implies that losses are caused by mainshocks while fore- and aftershocks yield no damage. However, recent earthquakes, such as the 2011 Tohoku, 2011 Christchurch, and 2016 Amatrice, have shown that aftershocks can also cause extensive damage, thus highlighting the need to study the effects earthquake sequences in loss model evaluation and cat bond pricing (Cousins et al., 2012; Stewart et al., 2018; Swiss Re, 2019). There are very few studies focusing on including these effects within the seismic risk assessment framework for spatially distributed assets. The few available studies modeled these effects using fitted model for aftershocks representation (Shokrabadi & Burton, 2019; Shome & Williams, 2014) or by using Epidemic-Type Aftershock Sequence (ETAS) model (Ogata, 1998; Zhang et al., 2018). Recently, Papadopoulos and Bazzurro (2021) studied the effects of seismicity clustering and damage accumulation on spatially distributed assets in the Umbria region in central Italy. Additionally, they developed damage-dependent fragility functions for the typical Italian building stocks. All the above studies concluded that neglecting the earthquake sequencing effects may lead to extensive underestimation of losses.

Building upon these studies, this study makes an attempt to propagate the earthquake sequencing effect in the
<figure><img src="image_4.png" /><figcaption>**FIGURE 2** Five municipalities located in central Italy selected for case-study implementation: (i) Perugia (dark blue); (ii) L'Aquila (red); (iii) Terni (light blue); (iv) Norcia (salmon red); (v) Amatrice (green).</figcaption></figure>

financial model, more specifically in the aggregate loss model, which is one of the key component for pricing risk-based cat bonds. The article presents a time-dependent aggregate loss model that accounts for the seismicity clustering and damage accumulation due to earthquake sequences. Moreover, we investigate the influence of time-dependency on the loss estimation and cat bond price. The structure of the remaining part of the article is as follows: First, we provide an overview of different components involved in the risk-based cat bond pricing framework followed by a description of the proposed time-dependent aggregate loss model. The article then shows an implementation of the proposed model using a case study analysis on five Italian municipalities. Finally, the influence of time-dependency on loss estimation and cat bond pricing is studied using a sensitivity analysis.

## 2 | OVERVIEW OF RISK-BASED CAT BOND PRICING

Pricing of risk-based cat bonds involves developing peril-specific catastrophe models: hazard, exposure, vulnerability, and financial (Hofer et al., 2020; Mistry & Lombardi, 2022, 2023). The hazard model addresses the following four questions: (i) What is the location of potential events? (ii) How large or severe the events are? (iii) What is the occurrence frequency of these events? (iv) What is the hazard intensity of event in the affected region? The first three questions are answered by generating multiple stochastic earthquake catalogs, which include large number of simulated earthquakes events representing broad spectrum of plausible events. For the last question, the hazard intensity for each event in the affected region is estimated using Ground Motion Model (GMM) along with local site conditions and other seismic parameters (such as magnitude, source to site

distance, style of faulting, inter- and intraevent variability). The exposure model defines key attributes of all the exposed assets, which include location, type of construction material, number of storeys, number of dwellings, economic activity type, built-up area, and replacement cost. The vulnerability model explains the damage caused to the asset given an occurrence of event with specific intensity level. These can be modeled using fragility and consequence functions. The former defines the probability of exceeding a damage state conditioned on the intensity measure (IM), while the latter defines range of damage ratios corresponding to each damage state. The convolution of hazard, exposure, and vulnerability model provide the cumulative distribution function for loss, which is further used as an input in the financial model. The financial model is comprised of three subcomponents: (i) interest rate model, (ii) aggregate loss model, and (iii) payoff functions. The interest rate model defines the rate of interest to be paid to the investors over different levels of maturity time. The aggregate loss model describes the frequency of events and their severity in terms of losses over time. The payoff function explains about the triggering mechanism of payment and the amount to be paid to sponsor and investor in case of occurrence/nonoccurrence of trigger event during the period of contract.

### 2.1 | Aggregate loss model

The aggregate loss model is a compound process with two underlying stochastic processes: (i) catastrophe events process, $N(t)$; (ii) catastrophe severity process, $X_n$. The former process simulates the occurrence of events at time, $t$, while the latter describes the severity of events at time, $t$. This model is based on the following three assumptions (Ma & Ma, 2013):
<figure><img src="image_5.png" /><figcaption>FIGURE 3 Representation of the workflow implemented for investigating the time-dependency in risk-based cat bond pricing framework.</figcaption></figure>

1. A Poisson point process $N(t)$ ($t \in [0, T]$) with intensity parameter $\lambda$, is used to describe the occurrence frequency of potential catastrophic event of a specific magnitude. The time instants of potential catastrophic events is denoted as $0 \le t_1 \le \dots \le t_n \le \dots \le T$. Here, $T$ denotes the total length of the cat bond contract.

2. The trigger mechanism is defined as the aggregate loss ($L_t$) exceeding a loss threshold level ($D$), mathematically the trigger event is expressed as: $\tau = \inf\{t: L(t) \ge D\}$.

3. The severity of catastrophic events at time $t_n$ is modeled as a random variable $\{X_n\}_{n=1,\dots}$, which is independent and identically distributed. The cumulative distribution function is expressed as $F(x) = P(X_n < x)$.

Mathematically, aggregate loss at a specific time, $L(t)$, can be expressed as:

$$L(t) = \sum_{n=1}^{N(t)} X_n. \qquad (1)$$

The cumulative distribution function of aggregate loss is expressed as:

$$F(D, T) = \sum_{n=0}^{\infty} e^{-\lambda T} \frac{(\lambda T)^n}{n!} F^n(D), \quad (2)$$

where $\lambda$ is the intensity parameter; $F^n(D)$ is the $n^{th}$ convolution of $F(x) = P(X_n < x)$. The intensity parameter $\lambda$ is calibrated using the temporal occurrence of earthquakes within a specific region. For the catastrophe severity process, $X_n$, the cumulative distribution function $F(x) = P(X_n < x)$ is calibrated using the simulated loss data. So far in previous studies (Hofer et al., 2020; Mistry & Lombardi, 2022, 2023), $\lambda$ calibration was performed using a homogeneous Poisson process while $F(x)$ calibration was based on the simulated loss data computed from time-independent seismic loss assessment. The main limitation of such an approach is the inability to account for seismicity clustering (foreshocks, mainshocks, and aftershocks) effects and damage accumulation due to earthquake sequences. In the current study, we eliminate these limitations by developing a time-dependent aggregate loss model. The details of proposed model is explained in the following section.

### 2.1.1 | Time-dependent aggregate loss model

We propose a time-dependent aggregate loss model where the intensity parameter of catastrophe event process is modeled using ETAS approach (Ogata, 1998) while the catastrophe severity process is modeled from simulated loss data computed using time-dependent seismic loss assessment.
<figure><img src="image_6.png" /><figcaption>**FIGURE 5** Shear wave velocity in the upper 30 m, $V_{s,30}$ with a spatial resolution of 50 m × 50 m. $V_{s,30}$ values are extracted from Mori et al. (2020).</figcaption></figure>

<figure><img src="image_7.png" /></figure>

The proposed model has the ability to account for seismicity clustering and damage accumulation due to earthquake sequences. The cumulative distribution function of the proposed aggregate loss is expressed as:

$$F(D, T)_{ES} = \sum_{n=0}^{\infty} e^{-\lambda_t T} \frac{(\lambda_t T)^n}{n!} F^n(D)_{ES}, \quad (3)$$

where $\lambda_t$ is the intensity parameter at a specific time, which is modeled using ETAS model; $F^n(D)_{ES}$ denotes the nth convolution of cumulative distribution function of losses calibrated using simulated loss data obtained from time-dependent seismic loss assessment. The ETAS model is based on the concept of self-exciting point-process where each

independent event has the ability to trigger an event. In this study, we use the ETAS model with exponentially tapered Omori time kernel (Mizrahi et al., 2021; Nandan et al., 2021; Ogata, 1998) and the intensity parameter at time t can be expressed as:

$$\lambda_t = \mu(x, y) + \sum_{i:t_i<t} g(t - t_i, x - x_i, y - y_i, m_i), \quad (4)$$

where $\mu$ represents intensity of background seismicity at specific location x, y; $t_i$, $x_i$, $y_i$, $m_i$ denotes the time, longitude, latitude, and magnitude of the ith event. Detailed description on the model and generation of catalogs are provided in Mizrahi et al. (2021).
<figure><img src="image_8.png" /><figcaption>**FIGURE 6** Time-independent seismic hazard map expressed in terms of peak ground acceleration (PGA) with 10% exceedance probability in 50 years (475-year return period).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 7** Time-dependent seismic hazard map expressed in terms of peak ground acceleration (PGA) with 10% exceedance probability in 50 years (475-year return period).</figcaption></figure>

## 3 | CASE-STUDY IMPLEMENTATION: ITALIAN MUNICIPALITY

In the present study, we focus on five municipalities located in central Italy: Perugia, L’Aquila, Terni, Norcia and Amatrice (see Figure 2). This case study has a twofold objective: (i) to illustrate the implementation of proposed time-dependent aggregate loss model and (ii) to study the impact of time-dependency on two loss metrics (annual exceedance probability and average annual loss, *AAL*) and cat bond price. Figure 3 illustrates the workflow employed in this study. First, we perform two different types of seismic loss assessment: time-independent (TI) and time-dependent (TD). Thereafter, the outcomes of both assessments are used as an input for developing the time-independent and time-dependent aggregate loss model, respectively.

### 3.1 | Hazard model

The first step in the hazard model is the generation of stochastic earthquake catalogs representing sufficiently long stochastic event set (say 100,000 years) to achieve convergence in hazard and risk metrics. In this study, we develop two hazard models: time-independent (TI) and time-dependent (TD) where each model includes 100,000 realization of 1-year-long stochastic earthquake catalog representing 100,000 years of seismicity (Musson, 2000). For TI model, the occurrence of earthquake is modeled as a homogeneous Poisson process (Crowley & Bommer, 2006) along with the earthquake source model (see Figure 4) and recurrence parameters used from the European Seismic Hazard Model (ESHM20) project (Danciu et al., 2021). The events generated with this model are independent in time, which
<figure><img src="image_10.png" /><figcaption>**FIGURE 8** Amplification in the seismic hazard map due to time-dependency. TI and TD refers to time-independent and time-dependent models.</figcaption></figure>

means two consecutive events in the catalog are independent of each other. For TD model, we employ ETAS approach with exponentially tapered Omori time kernel (Ogata, 1998) to generate independent and triggered earthquake events, which represents the dependency between two consecutive events over a time period (Mizrahi et al., 2021). The ETAS parameters were adopted from Hernandez et al. (2023) where we also present the validation of these parameters. These catalogs are then used to perform probabilistic seismic hazard analysis (PSHA) and seismic loss estimation. The former provides spatial distribution of hazard IM (peak ground acceleration, PGA, or spectral acceleration at first natural period, SaT₁) for a specific return period while the latter computes the distribution of loss levels with occurrence frequency. To estimate the ground-motion intensity at a specific location, we use the latest GMM proposed by Lanzano et al. (2019), which is calibrated using the updated Italian earthquake data set including the recent 2016–2017 central Italy earthquake sequence. Additionally, we account for inter- and intraevent uncertainty in the GMM by sampling εinter and εintra from standard normal distribution. The local site effect is taken into consideration by using the shear wave velocity in the upper 30 m, Vs,30 map for Italy (Mori et al., 2020) (see Figure 5). It is worth noting that spatial correlation is neglected in this study.

Figures 6 and 7 present the seismic hazard map for TI and TD models, respectively. These maps are expressed in terms of 10% probability of exceeding a PGA level in 50 years of time period. Figure 8 shows the amplification in seismic hazard due to the effect of time-dependency. Among the five municipalities, L’Aquila has the highest seismic hazard with PGA value of 0.5 g, which amplifies by the factor of 1.4 (0.7 g) in case of TD model. In Perugia, wider range of PGA levels are observed starting from 0.2 g (TI model) and 0.35 g (TD model), which goes as high as 0.4 g (TI model) and 0.6 g (TD model). Highest level of amplification (1.75) due to time-dependency is observed in eastern parts of Perugia while southwestern part observes lower amplification of 1.35.

In Amatrice and Norcia, for TI-model, PGA values range between 0.25 g and 0.4 g where higher seismic hazard is concentrated at the borders near the northern parts of Amatrice and southern parts of Norcia. These values are amplified with an average of 1.6 in case of TD model. In comparison with these municipalities, Terni exhibits lowest level of seismic hazard that ranges from 0.2 g to 0.35 g (TI model) and 0.35 g to 0.5 g (TD model). From the above observation, it can be concluded that neglecting the effect of time-dependency may lead to an average underestimation of 50% in the seismic hazard, which will lead to lower levels of loss estimates.

## 3.2 | Exposure model

We build a high-resolution grid-based exposure model by using spatial deaggregation approach (Dabbeek et al., 2021; Yepes-Estrada et al., 2017) where municipality level asset attribute is distributed onto a 100 m × 100 m (approximately, 3 arc-seconds) weighted grid. Figure 9 illustrates various steps involved in developing grid-based exposure model. The asset attribute data are extracted from European exposure data set repository (Crowley et al., 2020), which is derived from the latest Italian national census database (ISTAT, 2011). For the generation of weighted grid we used a 100-m resolution Global Human Settlement built-up layer 2022 (Pesaresi & Politis, 2022; Schiavina et al., 2022), which provides the spatial distribution of built-up area for a specific region (see Figure 10). It is worthy to note that only residential buildings were considered for the exposure model.

## 3.3 | Vulnerability model

Vulnerability model includes two subcomponents: (i) fragility function and (ii) consequence loss function. The former provides the probability of exceeding a specific damage
Gridded exposure modelling

<figure><img src="image_11.png" /><figcaption>**FIGURE 9** Graphical illustration of developing grid-based exposure model.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 10** Built-up area map for the selected municipality extracted from global human settlement layer (Pesaresi & Politis, 2022; Schiavina et al., 2022).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**FIGURE 11** Proportion of assets in each damage states: (A) TI model; (B) TD model; (C) relative difference in the percentage of assets between TD and TI models; M and RC denotes masonry (hashed bars) and reinforced concrete (solid bars).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**FIGURE 12** Annual exceedance probability of losses for TI model (solid line) and TD model (dashed line): (A) large municipalities: Perugia, L'Aquila, and Terni; (B) small municipality: Norcia and Amatrice.</figcaption></figure>

state conditioned on the hazard IM while the latter provides different levels of damage ratios for a specific damage state. In this study, we use suits of fragility and consequence function proposed in Papadopoulos and Bazzurro (2021).

For the time-independent seismic loss estimation, we use standard fragility functions, which are built assuming the initial damage state of the building as undamaged. On the other hand, for time-dependent seismic loss estimation, we
<table><caption>TABLE 1 Fragility curves corresponding to different vulnerability classes adopted from Papadopoulos and Bazzurro (2021).</caption><thead><tr><th>Material</th><th>Storeys</th><th>Vulnerability class</th><th>Description</th></tr></thead><tbody><tr><td rowspan="2">Masonry</td><td>1–2</td><td>URM_LR</td><td>Low-rise brick masonry</td></tr><tr><td>≥ 3</td><td>URM_MR</td><td>High-rise brick masonry</td></tr><tr><td rowspan="2">Reinforced concrete</td><td>1–3</td><td>RCINF_LR</td><td>Low-rise reinforce concrete</td></tr><tr><td>≥4</td><td>RCINF_MR</td><td>High-rise reinforce concrete</td></tr></tbody></table>

use damage-dependent fragility functions, which accounts for damage accumulation. Table 1 shows the fragility curves used for different vulnerability classes.

## 3.4 | Loss estimation

The loss estimation involves computing annual exceedance probability of losses and AAL. The former provides the probability of exceeding a loss value for a given return period. In simple terms it refers to per event occurrence which is typically known as occurrence exceedance probability (OEP). On the other hand, the AAL explains the occurrence of mean loss in any given year. These loss metrics are estimated using event loss tables, which is built from the convolution of hazard, exposure, and vulnerability model. In this study, we develop two event loss table each for time-independent (TI) and time-dependent (TD) loss estimation.

For the TI model, we use Monte Carlo simulation approach where the losses for ith event in TI stochastic earthquake catalog are calculated based on the assumption that initial damage state of the building is undamaged (Crowley & Bommer, 2006). Following are the steps involved in time-independent loss estimation:

(a) For ith event in jth TI stochastic earthquake catalog, calculate asset-specific IM using the GMM.

(b) Estimate the damage state of the asset by mapping the IM value (estimated in step 1) on the standard fragility function. Thereafter, damage ratio is estimated by sampling a random number between the minimum and maximum bounds of damage ratio for specific damage state.

(c) The loss is then calculated as a product of damage ratio, floor area (in sqm.) and replacement cost (per sqm.).

(d) Repeat steps 1–3 for all the events in jth catalog.

On the other hand, for TD model, a Markov chain Monte Carlo approach is used. In this case, the initial damage state of the asset is updated depending on the damage state of previous event (i – 1) (Papadopoulos & Bazzurro, 2021; Trevlopoulos & Guéguen, 2016; Trevlopoulos et al., 2020). The main advantage of this approach is the ability to account for damage accumulation over time due to earthquake sequences. Following are the steps involved in time-dependent loss estimation:

1. For first event (i = 1) in jth TD stochastic earthquake catalog, calculate asset-specific IM using the GMM.

2. Set the initial damage state $DS_{ini,(j,i)}$ and total damage ratio $DR_{total,(j,1)}$ as undamaged and 0, respectively.

3. Estimate current damage state $DS_{curr,(j,i)}$ and current damage ratio $DR_{curr,(j,i)}$ by considering the IM on the damage-dependent fragility function followed by sampling a random number between the minimum and maximum bounds of damage ratio corresponding to $DS_{curr,(j,i)}$. It is worth noting that event-specific damage ratio $DR_{event,(j,i)}$ for the first event in each catalog will be equal to $DR_{curr,(j,i)}$.

4. Now move on to the next event (i + 1) in jth stochastic earthquake catalog and repeat step 1.

5. Update the initial damage state and damage ratio as follows:

$$DS_{ini,(j,i+1)} = DS_{curr,(j,i)}, \quad (5)$$

$$DR_{total,(j,i+1)} = DR_{total,(j,i)} + DR_{event,(j,i)}. \quad (6)$$

6. Estimate the $DS_{curr,(j,i+1)}$ and $DR_{curr,(j,i+1)}$ similar to step 3 but in this case use the damage-dependent fragility function that represents the initial damage state of asset as $DS_{ini,(j,i+1)}$. The $DR_{event,(j,i+1)}$ is then estimated as follows:

$$DR_{event,(j,i+1)} = \begin{cases} DR_{new} - DR_{total,(j,i)} & ; DR_{total} \le DR_{curr} \\ 0 & ; DR_{total} > DR_{curr} \end{cases} \quad (7)$$

7. Finally, loss for each event are calculated as the product of event-specific damage ratio $DR_{event}$, floor area (in sqm.) and replacement cost (per sqm.).

8. Repeat steps 3–7 for all the events in the jth catalog.

### 3.4.1 | Influence of time-dependency on loss estimates

Figure 11A and B highlight the percentage of total assets in each damage states for TI and TD models, respectively. Figure 11C shows the relative difference in percentage of total assets between these models. It is worth noting that these results refer to the 475-year return period, which is estimated using the 100,000 years of simulation. In the case of TD model, a reduction of up to 20% is observed from the first two damage states (DS0 and DS1) while an increment of up to 38% is observed for the final damage state (DS4). Such a trend is observed because the TD model accounts for the damage accumulation over time due to earthquake sequences, which implicates that majority of the assets have moved from lower to higher damage states.

Figure 12A and B show the annual exceedance probability for large (Perugia, L’Aquila, and Terni) and small (Norcia and Amatrice) municipalities considered in the present case study, respectively. Each point on this curve represents the probability of exceeding a loss value in a given
<figure><img src="image_15.png" /><figcaption>**FIGURE 14** TD model: spatial distribution of losses (in €million) for 475-year return period.</figcaption></figure>

<figure><img src="image_16.png" /></figure>

return period. It is observed that the TD model for all the municipalities exhibits higher loss levels as compared to the TI model. The influence of time-dependency on these loss estimates is further investigated by comparing the spatial distribution of loss value corresponding to 475-year return period (10% exceedance probability in 50 years) for both TI and TD models (see Figures 13 and 14, respectively). Figure 15 presents the amplification due to time-dependency on the loss value corresponding to 475-year return period. For Perugia and Terni, densely populated area experiences amplification ranging between 1.3 and 2.0 while lower amplification ranging between 0.5 and 1.5 are observed for other municipalities (i.e., L'Aquila, Norcia, and Amatrice). Moreover, higher amplification range (3.0–4.5) is observed in the southwestern parts of Perugia and Terni with few sites experi-

encing amplification as high as 7. For both TI and TD models, a common observation can be made that higher loss levels are concentrated in the densely populated region for all the municipalities.

Figure 16A shows the AAL for TI and TD models, which is computed as the ratio of sum of losses from each catalog to the total number of catalog years in the event loss table. Among the five municipalities, for TI model, Perugia experiences highest AAL estimate (€15.6 mn) followed by L'Aquila (€11.3 mn), Terni (€8.78 mn), Norcia (€0.81 mn), and Amatrice (€0.55 mn). In case of TD model, an average amplification of 3.64 is observed in the AAL estimates for all the municipality. Figure 16B presents the two sets of cumulative distribution function for losses for TI (solid lines) and TD models (dashed lines), which are further used as an
<figure><img src="image_17.png" /><figcaption>**FIGURE 15** Amplification in the loss values due to time-dependency. TI and TD refers to time-independent and time-dependent models.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**FIGURE 16** (A) Average annual losses for TI model (light gray) and TD model (dark green). AF denotes amplification factor. (B) cumulative distribution function for losses: TI model (solid line) and TD model (dashed line).</figcaption></figure>

input for respective aggregate loss model. These functions are obtained by fitting a lognormal distribution to the simulated loss data. The estimated mean and standard deviation of lognormal distribution are presented in the Appendix (see Appendix A).

**3.5 | Pricing risk-based cat bonds**

The price of cat bonds are estimated by combining the interest rate model, aggregate loss model, and payoff function. In this study, we use the Cox–Ingersoll–Ross (CIR) model (Cox et al., 1985) as the interest rate model along with the zero-coupon payoff function. Detailed pricing equations are provided in the Appendix (see Appendix B). We calibrated the CIR model parameters ($\theta = 3.03\%$, $\sigma = 5.59\%$, $k = 0.0533$, $\lambda_r = -0.01$, $r_0 = 3.03\%$) by using 3-month

maturity US monthly treasury bill data ranging between 1994 and 2022. For both TI and TD models, cat bonds are priced at face value of €1 at time $t = 0$ for different levels of time to maturity, $T \in [0.25 \text{ years}, 3 \text{ years}]$ and loss thresholds, $D \in [€0.01 \text{ mn}, €5000 \text{ mn}]$. In case of TI aggregate loss model, a single intensity parameter ($\lambda = 0.5$) for Equation (2) is estimated as the mean intensity parameter of all the seismic source zone from the hazard model. On the other hand, for TD aggregate loss model, the intensity parameter for each level of time to maturity $T$ is calculated using Equation (4). The summary of steps involved in estimating intensity parameter for TD model are as follows:

1. Generate 100,000 realization of 3-year-long TD stochastic earthquake catalog using the ETAS model, which represents 100,000 years of seismicity for three consecutive years.
<figure><img src="image_19.png" /><figcaption>**FIGURE 17** Zero-coupon catastrophe bond price for TI and TD models. The contour lines represent different levels of time to maturity. TI and TD denote time-independent and time-dependent models.</figcaption></figure>

2. For each realization split the 3-year-long stochastic catalog into 36 subcatalogs representing each month (3 years = 36 months).

3. Compute the total intensity at *i*th time to maturity *T*ᵢ by adding the background intensity rate with the summa-

tion of aftershock intensity rate from all the events that occurred until time *T*ᵢ.

4. Finally, repeat steps 2–3 for all the 100,000 realization and take the mean of total intensity value for each month.
<figure><img src="image_20.png" /><figcaption>**FIGURE 18** Relative difference in cat bond price between TI and TD model for three levels of time to maturity: *T* = 1, 2, 3 years. TI and TD denotes time-independent and time-dependent models.</figcaption></figure>

**3.6 | Influence of time-dependency on cat bond price**

Figure 17 presents the contour plot of cat bond price corresponding to given levels of loss threshold and time to maturity for TI and TD models. Two common observations can be made from both models: (i) at a given level of time to maturity, the price of cat bond is directly proportional to the loss threshold levels; (ii) at a given loss threshold level, the price of cat bond is indirectly proportional to the time to maturity levels. For investigating the influence of time-dependency on cat bond prices, we compute the relative difference between the TI and TD models (see Figure 18). Highest level of difference is observed at lower loss threshold levels and higher time to maturity. For example, if we compare the difference at a single time to maturity (say *T* = 3 years), then these differences increase as we move from low to high loss threshold (5% to 45%). On the other hand, at a single loss threshold level (say *D* = €50 mn), the difference increases from low to high time to maturity (25% to 45%). Additionally, when we compare municipalities, higher differences are observed in case of large municipalities. This implicates that large municipalities have higher level of influence due to time-dependency as compared to the small one. From the above observation we can conclude that neglecting time-dependency in the cat bond pricing leads to overestimation of up to 45% at lower loss thresholds and up to 5% for higher threshold.

**4 | CONCLUSION**

This study presents a new time-dependent aggregate loss model (TD model). In the proposed approach, the occurrence of catastrophe events is modeled using an epidemic-type aftershock sequencing approach while the severity of the events is obtained from simulated loss data computed using

the Markov chain Monte Carlo algorithm. The implementation of the proposed method is presented as a case study analysis on five Italian municipalities; results are then compared with those obtained from a time-independent model. The comparison shows that the use of the TD model has a significant influence on the distribution of the hazard, computed losses, and cat bond pricing. More specifically, results obtained from the proposed TD model result in an amplification between 1.3 and 1.75 in the seismic hazard map corresponding to 475-year return period. In terms of damage states, a higher proportion of assets moves from low damage state to the higher damage state. Losses associated with 475-year return period show an amplification up to 2.0, while an average of 3.4 amplification is observed for AAL. The time-dependency causes reduction (up to 45%) in the cat bond price at lower threshold levels and longer times to maturity. Future research is required to investigate the impact of interest rates, spatial correlation between hazard intensity, and damage levels. Finally, the proposed model can be adapted to account for losses caused by cascading events from multiple perils.

**AUTHOR CONTRIBUTIONS**

H.K.M. initially drafted the paper, developed the algorithm and models, performed the analysis, and interpreted the results. A.H. and P.G. provided guidance with the conceptual understanding of time-dependent model (ETAS) and result interpretations. A.H. estimated the ETAS parameters for the study region. D.L. provided guidance for the development of time-independent (Poisson) modeling algorithm and result interpretation. All the authors reviewed and approved the final draft.

**ACKNOWLEDGMENTS**

The first author wishes to acknowledge the funding received from *MACE Prize Doctoral Scholarship* and *Turing Scheme: PGR research exchange* in support of his doctoral studies.
All the simulation for this work was carried out using High Performance Computing cluster—Computational Shared Facility 3 (CSF3) at The University of Manchester. A.H. is funded by the URBASIS-EU project (H2020-MSCA-ITN-2018), Grant 813,137. P.G. would like to thank LabEx OSUG@2020 (Investissements d'avenir-ANR10LABX56) and the AXA Research Fund supporting the project New Probabilistic Assessment of Seismic Hazard, Losses and Risks in Strong Seismic Prone Regions.

## CONFLICT OF INTEREST STATEMENT

The authors declare no potential conflict of interests.

## ORCID

Harsh K. Mistry https://orcid.org/0000-0003-0593-0605

## REFERENCES

Allen, R. M., & Melgar, D. (2019). Earthquake early warning: Advances, scientific challenges, and societal needs. *Annual Review of Earth and Planetary Sciences*, *47*, 361–388.

Bayliss, C., Guidotti, R., Estrada-Moreno, A., Franco, G., & Juan, A. A. (2020). A biased-randomized algorithm for optimizing efficiency in parametric earthquake (RE) insurance solutions. *Computers & Operations Research*, *123*, 105033.

Burnecki, K., & Kukla, G. (2003). Pricing of zero-coupon and coupon cat bonds. *Applicationes Mathematicae*, *30*, 315–324.

Cardenas, V., Hochrainer, S., Mechler, R., Pflug, G., & Linnerooth-Bayer, J. (2007). Sovereign financial disaster risk management: The case of Mexico. *Environmental Hazards*, *7*(1), 40–53.

Cochran, E. S., Kohler, M. D., Given, D. D., Guiwits, S., Andrews, J., Meier, M.-A., Ahmad, M., Henson, I., Hartog, R., & Smith, D. (2018). Earthquake early warning shakelimit system: Testing and certification platform. *Seismological Research Letters*, *89*(1), 108–117.

Cousins, W., King, A., & Kanga, M. (2012). Accumulated losses from sequences of earthquakes: implications for risk modeling. In *Proceedings of the 15th world conference on earthquake engineering* (Vol. 24, p. 28). WCEE.

Cox, J. C., Ingersoll Jr, J. E., & Ross, S. A. (1985). A theory of the term structure of interest rates. *Econometrica*, *53*(2), 385–407.

Crowley, H., & Bommer, J. J. (2006). Modelling seismic hazard in earthquake loss models with spatially distributed exposure. *Bulletin of Earthquake Engineering*, *4*(3), 249–273.

Crowley, H., Despotaki, V., Rodrigues, D., Silva, V., Costa, C., Toma-Danila, D., Riga, E., Karatzetzou, A., Fotopoulou, S., Sousa, L., Ozcebe, S., Gamba, P., DabbEEK, J., Romão, X., Pereira, N., Castro, J. M., Daniell, J., Veliu, E., Bilgin, H., ... Contributors, E. (2020). *European exposure model data repository (version 1.0) [data set]*. http://doi.org/10.5281/zenodo.4062044

DabbEEK, J., Crowley, H., Silva, V., Weatherill, G., Paul, N., & Nievas, C. I. (2021). Impact of exposure spatial resolution on seismic loss estimates in regional portfolios. *Bulletin of Earthquake Engineering*, *19*(14), 5819–5841.

Danciu, L., Nandan, S., Reyes, C., Basili, R., Weatherill, G., Beauval, C., Rovida, A., Vilanova, S., Sesetyan, K., Bard, P., Cotton, F., & Wiemer, S. (2021). *The 2020 update of the European Seismic Hazard Model: Model overview*. EFEHR Tech Rep, 1, v1.

Franco, G. (2010). Minimization of trigger error in cat-in-a-box parametric earthquake catastrophe bonds with an application to Costa Rica. *Earthquake Spectra*, *26*(4), 983–998.

Härdle, W. K., & Cabrera, B. L. (2010). Calibrating cat bonds for Mexican earthquakes. *Journal of Risk and Insurance*, *77*(3), 625–650.

Hernandez, A., Mistry, H. K., Gueguen, P., & Trevlopoulos, K. (2023). Development of a time-dependent vulnerability model for aftershock sequences. (in prep).

Hofer, L., Gardoni, P., & Zanini, M. A. (2019). Risk-based cat bond pricing considering parameter uncertainties. *Sustainable and Resilient Infrastructure*, *6*(5), 315–329.

Hofer, L., Zanini, M. A., & Gardoni, P. (2020). Risk-based catastrophe bond design for a spatially distributed portfolio. *Structural Safety*, *83*, 101908.

ISTAT. (2011). https://www.istat.it/en/

Kajwang, B. (2020). Factors influencing the low penetration of insurance services in Africa. *Journal of Developing Country Studies*, *5*(1), 25–35.

Kelly, M., Bowen, S. G., & McGillivray, R. G. (2020). The earthquake insurance protection gap: Tale of two countries. *Journal of Insurance Regulation*, *20*(11), 1–37.

Lanzano, G., Luzi, L., Pacor, F., Felicetta, C., Puglia, R., Sgobba, S., & D'Amico, M. (2019). A revised ground-motion prediction model for shallow crustal earthquakes in Italy revised ground-motion prediction model for shallow crustal earthquakes in Italy. *Bulletin of the Seismological Society of America*, *109*(2), 525–540.

Ma, Z.-G., & Ma, C.-Q. (2013). Pricing catastrophe risk bonds: A mixed approximation method. *Insurance: Mathematics and Economics*, *52*(2), 243–254.

Mistry, H. K., & Lombardi, D. (2022). Pricing risk-based catastrophe bonds for earthquakes at an urban scale. *Scientific Reports*, *12*(1), 1–12.

Mistry, H. K., & Lombardi, D. (2023). A stochastic exposure model for seismic risk assessment and pricing of catastrophe bonds. *Natural Hazards*, *117*, 803–829.

Mizrahi, L., Nandan, S., & Wiemer, S. (2021). The effect of declustering on the size distribution of mainshocks. *Seismological Society of America*, *92*(4), 2333–2342.

Mori, F., Mendicelli, A., Moscatelli, M., Romagnoli, G., Peronace, E., & Naso, G. (2020). A new vs30 map for Italy based on the seismic microzonation dataset. *Engineering Geology*, *275*, 105745.

Musson, R. M. (2000). The use of monte carlo simulations for seismic hazard assessment in the uk. *Annali Di Geofisica*, *43*(1), 1–9.

Nandan, S., Kamer, Y., Ouillon, G., Hiemer, S., & Sornette, D. (2021). Global models for short-term earthquake forecasting and predictive skill assessment. *The European Physical Journal Special Topics*, *230*(1), 425–449.

Ogata, Y. (1998). Space-time point-process models for earthquake occurrences. *Annals of the Institute of Statistical Mathematics*, *50*(2), 379–402.

Papadopoulos, A. N., & Bazzurro, P. (2021). Exploring probabilistic seismic risk assessment accounting for seismicity clustering and damage accumulation: Part II. risk analysis. *Earthquake Spectra*, *37*(1), 386–408.

Pesaresi, M., & Politis, P. (2022). *GHS built-up surface grid, derived from sentinel2 composite and landsat, multitemporal (1975-2030)*. European Commission, Joint Research Centre (JRC).

Pothon, A., Gueguen, P., Buisine, S., & Bard, P.-Y. (2019). California earthquake insurance unpopularity: The issue is the price, not the risk perception. *Natural Hazards and Earth System Sciences*, *19*(8), 1909–1924.

Schiavina, M., Melchiorri, M., Pesaresi, M., Politis, P., Freire, S., Maffenini, L., Florio, P., Ehrlich, D., Goch, K., Tommasi, P., & Kemper, T. (2022). *GHSF data package 2022*. European Commission, Joint Research Centre (JRC).

Shao, J., Pantelous, A. A., Ayyub, B. M., Chan, S., & Nadarajah, S. (2017). Nuclear catastrophe risk bonds in a Markov-dependent environment. *ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part A: Civil Engineering*, *3*(4), 04017018.

Shao, J., Papaioannou, A., & Pantelous, A. (2016). Pricing and simulating cat bonds in a Markov-dependent environment. *Applied Mathematics and Computation*, *309*, 68–84.

Shokrabadi, M., & Burton, H. V. (2019). Regional short-term and long-term risk and loss assessment under sequential seismic events. *Engineering Structures*, *185*, 366–376.

Shome, N., & Williams, C. (2014). Aftershock risk in Japan following Tohoku earthquake. In *Proceedings of the 10th US national conference on earthquake engineering, Alaska, Earthquake Engineering Research Institute* (pp. 21–25).
Stewart, J. P., Zimmaro, P., Lanzo, G., Mazzoni, S., Ausilio, E., Aversa, S., Bozzoni, F., Cairo, R., Capatti, M. C., Castiglia, M., Chiabrando, F., Chiaradonna, A., d'ONOFRIO, A., Dashti, S., De Risi, R., de Silva, F., Pasqua, F. D., Dezi, F., Di Domenica, A., ... Tropeano, G. (2018). Reconciliation of 2016 central Italy earthquake sequence. *Earthquake Spectra*, 34(4), 1547–1555.

Swiss Re, S. (2019). *Sigma: L'aquila 10 years on*. https://www.swissre.com

Swiss Re, S. (2023). *Sigma explorer: Catastrophe losses*. https://www.sigma-explorer.com

Trevlopoulos, K., & Guéguen, P. (2016). Period elongation-based framework for operative assessment of the variation of seismic vulnerability of reinforced concrete buildings during aftershock sequences. *Soil Dynamics and Earthquake Engineering*, 84, 224–237.

Trevlopoulos, K., Guéguen, P., Helmstetter, A., & Cotton, F. (2020). Earthquake risk in reinforced concrete buildings during aftershock sequences based on period elongation and operational earthquake forecasting. *Structural Safety*, 84, 101922.

Yepes-Estrada, C., Silva, V., Valcárcel, J., Acevedo, A. B., Tarque, N., Hube, M. A., Coronel, G., & Santa María, H. (2017). Modeling the residential building inventory in South America for seismic risk assessment. *Earthquake Spectra*, 33(1), 299–322.

Zhang, L., Werner, M. J., & Goda, K. (2018). Spatiotemporal seismic hazard and risk assessment of aftershocks of m 9 megathrust earthquakes. *Bulletin of the Seismological Society of America*, 108(6), 3313–3335.

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Mistry, H. K., Hernandez, A., Guéguen, P., & Lombardi, D. (2024). Effect of earthquake sequences on risk-based catastrophe bond pricing. *Risk Analysis*, 1–16.
https://doi.org/10.1111/risa.14288

## APPENDIX A: LOSS FUNCTIONS
### PARAMETERS (MEAN AND STANDARD DEVIATION)

Table A.1 presents all the parameters estimated by fitting lognormal distribution to the simulated loss data in Section 3.4.1.

## APPENDIX B: PRICING FORMULATION FOR CAT BONDS

The CIR model is a type of interest rate model, which computes the rate of interest to be paid to the bond holder at a specific time to maturity (Cox et al., 1985). Mathematically, assuming that the spot interest rate follows a mean-reverting square-root process, the rate of interest at time *t* can be expressed as:

$$dr(t) = k[\theta - r(t)]dt + \sigma\sqrt{r(t)}dW(t), \quad (B.1)$$

where *k* denotes the mean reverting force, *θ* is the long-term mean, *σ* refers to instantaneous volatility, *W(t)* is the

<table><caption>**TABLE A.1** Mean and standard deviation for lognormally fitted loss functions for time-independent (TI model) and time-dependent models (TD model).</caption><thead><tr><th>Municipality</th><th>Model</th><th>Mean</th><th>Std</th></tr></thead><tbody><tr><td rowspan="2">Perugia</td><td>TI-model</td><td>1.64</td><td>3.45</td></tr><tr><td>TD-model</td><td>3.54</td><td>3.46</td></tr><tr><td rowspan="2">L'Aquila</td><td>TI-model</td><td>1.3</td><td>3.4</td></tr><tr><td>TD-model</td><td>3.2</td><td>3.4</td></tr><tr><td rowspan="2">Terni</td><td>TI-model</td><td>1.44</td><td>3.31</td></tr><tr><td>TD-model</td><td>2.81</td><td>3.24</td></tr><tr><td rowspan="2">Norcia</td><td>TI-model</td><td>0.44</td><td>2.95</td></tr><tr><td>TD-model</td><td>0.88</td><td>2.92</td></tr><tr><td rowspan="2">Amarice</td><td>TI-model</td><td>0.31</td><td>2.95</td></tr><tr><td>TD-model</td><td>0.64</td><td>2.88</td></tr></tbody></table>

standard Weiner process where *t* is in range 0–T. The main advantage of this interest rate model is its ability to provide nonnegative interest rates, which is ensured using the following condition:

$$2k\theta > \sigma^2. \quad (B.2)$$

The pure discount T-bond at time *t* can be computed using the above conditions as follows:

$$B_{CIR}(t, T) = A(t, T)e^{-B(t,T)r(t)}, \quad (B.3)$$

$$A(t, T) = \left[ \frac{2\gamma e^{(k+\lambda_r+h)(T-t)/2}}{2\gamma + (k + \lambda_r + \gamma)(e^{(T-t)\gamma} - 1)} \right]^{2k\theta/\sigma^2}, \quad (B.4)$$

$$B(t, T) = \left[ \frac{2(e^{(T-t)\gamma} - 1)}{2\gamma + (k + \lambda_r + \gamma)(e^{(T-t)\gamma} - 1)} \right], \quad (B.5)$$

$$\gamma = \sqrt{(k + \lambda_r)^2 + 2\sigma^2}. \quad (B.6)$$

The zero-coupon payoff function can be represented as

$$P_{ZC} = \begin{cases} Z, & L(T) \le D \\ \eta Z, & L(T) > D, \end{cases} \quad (B.7)$$

where *Z* denotes the face value to be paid to bond holder in case of nontriggering of bond contract until the maturity time, *T*. The trigger event is defined as any event with aggregate loss exceeding a predefined loss threshold level. Finally, the zero-coupon cat bond can be priced under the risk-neutral pricing measure *Q* by combining the interest rate model, aggregate loss model, and payoff function. Mathematically it can be represented as:

$$ \begin{aligned} V_{ZC}(t) &= E^Q \left( Ze^{-\int_{t}^{T} r_s ds} P_{ZC} \| F_t \right) \\ &= B_{CIR}(t, T) Z \{ [F(D, T) + \eta[1 - F(D, T)]] \} \end{aligned} \quad (B.8) $$