A stochastic exposure model for seismic risk assessment and pricing of catastrophe bonds

Risk-based catastrophe bonds require the estimation of losses from the convolution of hazard, exposure and vulnerability models. These models are affected by different uncertainties that arise from the definition of their input parameters. In this paper, we propose a stochastic approach to treat the uncertainty in the asset location and attributes of the exposure model. The proposed method uses the Monte Carlo sampling approach to generate a stochastic exposure database, where each asset location is generated randomly within the geometric bound of the administration, while the asset attributes (i.e. construction type and material, number of storey, building activity type and number of dwelling) are sampled from distributions built from census data. Finally, a sensitivity analysis is performed to investigate the influence of the spatial resolution of the exposure model on the average annual losses (AAL) and catastrophe bond prices. To this end, we implement four exposure models, with spatial resolution at the asset, municipality and province levels, on a study region comprising ten provinces in southern Italy. Compared to the proposed model, the exposure model where assets are relocated and aggregated at the geometric centroid of the municipality underestimates AAL by 8%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8\%$$\end{document}, while a higher difference (up to 20%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$20\%$$\end{document}) is observed for the exposure model where assets are relocated and aggregated at the geometric centroid of the province. We also consider an exposure model whose asset locations are extracted from publicly available building footprints. Yet, this latter model was incomplete for some provinces resulting in underestimation of AAL up to 90%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$90\%$$\end{document}. Differences in catastrophe bond prices obtained from the four exposure models are less evident, with the exposure model built based on the building footprints showing a difference up to 9.5%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$9.5\%$$\end{document}.


Introduction
Earthquakes have a long history of causing severe human and economic losses. Both human and economic losses can be reduced by constructing earthquake-resistant buildings as well as updating current construction practises. Additionally in case of post-disaster recovery, economic losses can be managed by implementing disaster risk management strategies that transfer part of the financial risk from the impacted communities to other financial entities, typically re/insurance firms and the financial market. Such measures will help in boosting the post-disaster recovery of the affected community. Insurance and catastrophe bonds are the main financial products used in disaster risk management strategies. Earthquake insurance policies are typically underwritten by private individuals, such as asset owners, who seek financial protection against losses caused directly, or indirectly, by earthquakes in exchange of a regular payment, the premium. Catastrophe bonds, cat bonds in short, are a type of insurance-linked security issued by a sponsor, typically a re/ insurance firm or a governmental entity, to raise money and transfer part of the risk to the financial market. Compared to other bond types, cat bonds have short maturity dates, between 3 and 5 years, and are attractive to investors, such as institutional investors and hedge funds, who are looking for competitive returns and investment diversification. The popularity of cat bonds has increased over the past decade, reaching a new record in 2021 with $ 12.8 billion of new issuance (source: Bloomberg). Despite the maturity of these financial products, the majority of economic losses caused by earthquakes remain largely uninsured. Figure 1 presents the costliest earthquakes to the insurance industry between the period 1980 and 2020. Except for the events in New Zealand, all other have low insurance Fig. 1 Ten most expensive earthquakes to the insurance industry worldwide between 1980 and 2020 (Statista 2022) 1 3 penetration. Some of the main factors of low insurance penetration includes affordability, low risk perception, design and structure of the cat bonds (Bevere and Grollimund 2012; Kelly et al. 2020). These factors are influenced by uncertainty and low confidence in the loss estimates, which tend to increase the price of insurance policies, making them less attractive to potential customers. Although several studies on earthquake insurance have been published in the last decade (Yucemen 2005; Petseti and Nektarios 2012; Goda and Yoshikawa 2012;Asprone et al. 2013;Yoshikawa and Goda 2014;Tian and Yao 2015;Lin 2018;Goda et al. 2020;Gkimprixis et al. 2021), research on cat bonds has been limited (Burnecki and Kukla 2003;Härdle and Cabrera 2010;Ma and Ma 2013;Shao et al. 2016Shao et al. , 2017Hofer et al. 2019Hofer et al. , 2020 and mainly focused on pricing cat bond based on worldwide historical loss data sets or simulated loss data computed at low resolution. One of the early works on cat bonds used a doubly stochastic Poisson process and 10 years of catastrophe loss data to price zero-coupon and coupon bonds (Burnecki and Kukla 2003). Härdle and Cabrera (2010) priced cat bonds for Mexico using a hybrid (modelledindex) type of trigger mechanism that combined losses and moment magnitude of the seismic events. Ma and Ma (2013) investigated the effect of uncertainty in the interest rate when pricing cat bonds. Shao et al. (2016) proposed a new formulation to price cat bonds based on different types of payoff functions. Shao et al. (2017) developed a nuclear catastrophe risk bond using a Markov-dependent environment and the Vasicek and Cox-Ingersoll-Ross (CIR) interest rate model. Hofer et al. (2019) studied the uncertainty associated with claim arrival and claim severity rates, and introduced a new formulation to price cat bonds for different levels of accepted risk. Hofer et al. (2020) proposed a method to price cat bonds at a national level using a simplified hazard model and asset-independent replacement costs. A recent work by Mistry and Lombardi (2022) focused on cat bond pricing for a small urban area located in southern Italy.
Similarly to the pricing of insurance policies, cat bonds are priced based on aggregate losses. Aggregate losses can be computed using either historical loss data sets (Shao et al. 2016;Ma and Ma 2013) or simulated loss data (Hofer et al. 2020;Mistry and Lombardi 2022). The latter approach is known as risk-based cat bond pricing and requires the formulation and convolution of hazard, exposure, vulnerability models. Each model is affected by two types of uncertainty: aleatory and epistemic. The aleatory uncertainty arises from the inherent randomness of earthquakes (e.g. frequency and severity of events). On the other hand, the epistemic uncertainty is related to lack of knowledge and simplifications made in the modelling process. Moreover, uncertainty associated with exposure model is mainly due to limited information available regarding asset attributes (e.g. their occupancy, construction type, number of dwelling and replacement cost). It is worth noting that the uncertainty in vulnerability component (record-to-record, buildingto-building and uncertainty in the damage criterion) can also impact the risk results by 5-15%, but including it in the present framework was out of the scope. Interested readers may refer to some of the recent study by Sousa et al. (2016Sousa et al. ( , 2018; Silva (2019) that presented various approaches to treat these uncertainty and their implications on different loss metrics. Risk-based cat bond pricing requires that these uncertainties are quantified and taken into account explicitly. The research on cat bonds carried out so far has been limited, and mostly focused on financial aspects, such as interest rate model, types of bond triggering mechanisms.
Previous studies on seismic risk assessment have focused on treating the uncertainty in the hazard Crowley et al. 2008;Silva 2018) and vulnerability models (Sousa et al. 2016(Sousa et al. , 2018Silva 2019), with only a few studies treating the uncertainty in the exposure model. Studies focusing on treating the uncertainty in the exposure model mainly involve using different spatial resolution of aggregated exposure along with mapping schemes to translate the descriptors from national census data to building classes (Bal et al. 2010;Yepes-Estrada et al. 2017). Kalakonas et al. (2020) treated the exposure uncertainty by adopting five different spatial resolution for asset disaggregation and two alternative models for distributions of construction material. Dabbeek et al. (2021) quantified the impact of using exposure models with different spatial resolution on the seismic losses and recommended using grid resolution ranging between 120 and 240 arc-seconds in order to keep the bias within 5%. Differently from previous studies, we propose a new approach to generate a stochastic exposure model that accounts for the uncertainty in the asset location and its attributes. According to the proposed method, the asset location is generated by Monte Carlo sampling within the administrative bounds, while the asset attributes such as construction type and material, number of storey, building activity type and number of dwelling, are sampled from distributions computed from census data. Such sampling is repeated for a large number of realizations to account for the above-mentioned uncertainties in the exposure model. Furthermore, we perform a sensitivity analysis to investigate the effect of different spatial resolutions in the exposure model on the average annual losses (AAL) and cat bond pricing.

Overview of current databases and models
There have been several international efforts to develop global and national exposure models. The first open-access global exposure model was developed by the U.S. Geological Survey's Prompt Assessment of Global Earthquakes for Response (PAGER) programme (Jaiswal et al. 2010). The model was constructed by harmonizing data sets at a national scale from different sources, including the UN statistics, UN Habitat's demographic and health survey, national housing censuses data. The building classification was based on the construction material, lateral force resisting system, type of occupancy and location (i.e. urban or rural). De Bono and Mora (2014) developed an open-access global exposure model with a 5 km spatial resolution, which was released in the Global Assessment Report (GEG-2013). The Global Earthquake Model (GEM) foundation introduced the Global Exposure Database (GED4GEM) (Gamba 2014) for spatial distribution of residential buildings and population at three different spatial resolutions: national, administrative boundaries 1 and 2. The latest effort, led by Crowley et al. (2020), introduced the European Exposure model for the spatial distribution of buildings (residential, commercial and industrial), population and replacement costs for 45 countries.
Most of the available exposure models cited earlier provide sufficiently accurate information for risk assessment at a national and global scales. However, they tend to lack the resolution required to perform risk assessment at the urban scale with a high confidence level. The main challenge in developing a high-resolution exposure model lies in the limited information available at the asset level. For large urban areas, exposures data can be retrieved from open-source and proprietary databases, such as OpenStreetMap (OSM) (OpenStreetMap 2017) and Google Earth (Gorelick et al. 2017). However, such data are lacking, or incomplete, in small urban areas and rural regions. As a result, several simplifying methods have been proposed to model exposure data when information is limited, or not available. The available methods propose relocation and aggregation of assets at the centroid of administrative areas (Bazzurro and Park 2007;Bal et al. 2010), or adoption of weighted grids determined from population density data (Scheingraber and Käser 2019;Dabbeek et al. 2021), such as those provided by the Global Health Settlement Layer (Pesaresi et al. 2015). Bazzurro and Park (2007) proposed a method for relocation and aggregation of assets based on zip codes; yet, it was found that such an approach introduced artificial asset-to-asset correlation and led to an underestimation of losses arising from frequent events and overestimation of losses from extreme rare events. Bal et al. (2010) investigated the role played by the exposure model with four resolutions from coarse to fine (i.e. district, postcode, sub-district and geo-cell level). The study concluded that the exposure model with coarse resolution at district level provided unreliable estimates for losses. Scheingraber and Käser (2019) built an exposure model using weighted grid points within administrative zones, with weights determined based on the global population database LandScan (Dobson et al. 2000), digital elevation models and land cover data. Additionally, they studied the impact of location uncertainty on the probable maximum loss and found that neglecting location uncertainty led to significant inaccuracy in the results, especially for smaller portfolios. Dabbeek et al. (2021) performed a sensitivity analysis to study the effects of exposure spatial resolution on the seismic losses at regional level in 35 countries. In total, 18 test cases were analysed using different exposure modelling approaches (i.e. geometric centroid, density-weighted centroid and location of maximum density) and resolutions (i.e. admin levels and grid resolutions). The study revealed that the methods adopting relocation and aggregation of assets at the geometric centroid of the administrative boundaries led to bias of up to 27% in the average annual losses.

Proposed stochastic exposure model
We propose a method to generate a stochastic exposure model that uses the Monte Carlo sampling method to account for the asset location and its attributes-related uncertainty in the exposure model. Figure 2 illustrates the steps required to generate stochastic exposure model. The asset attributes, such as building material (BM), building storey (BS), building activity (BA) and number of dwelling in a building (DW), are modelled as random variables with known probability distributions. These distributions can be determined from available databases, such as national census data.
The proposed exposure model can be built in four steps (see also Fig. 2): (1) for each i th asset, the geographical coordinates (latitude and longitude) are randomly generated from a uniform distribution within the administrative bound.
(2) once the asset location is determined, a set of numbers between 0 and 1 are randomly generated from a standard uniform distributions representing the different asset attributes, including: construction material bm rand,i , number of storeys bs rand,i , type of activity ba rand,i , number of dwellings dw rand,i . Additionally, a relationship between construction material and number of storeys should be established. For example, if the sampled construction material for i th asset is masonry, then the maximum number of storey that can be sampled for this asset should be limited by the maximum number of storey a masonry building can have within its respective administrative bounds. It is worth noting that the attribute distributions and limits can be determined from census 1 3 data. The total building area can be computed as the product of the number of dwellings and the average floor area per dwelling. (3) depending on the type of activity, a replacement cost rc rand,i is sampled from a distribution, typically a uniform distribution bounded by a lower and upper replacement costs for a specific building activity. (4) the process is iterated until the total number of assets set out in the census data is reached. Finally, the data are compiled in a database which constitutes the stochastic exposure model.
In total N number of stochastic exposure models are generated, where N refers to the total number of earthquake events in the synthetic catalogue. The main advantage of the proposed model over the available exposure models discussed in the previous section is its ability to incorporate the uncertainty in the asset location and attributes, especially in the replacement costs, which can vary greatly from asset to asset even within a small urban area. The detailed implementation of the proposed model is explained in Sect. 4.4. There are two main limitations of the proposed method. First limitation is the assumption of sampling asset location randomly from uniformly distribution within the administrative bounds. This assumption is insufficient to exclude the locations that will not be urbanized, such as water bodies, reserves and other protected regions. For the future studies such limitation can be eliminated by sampling a location based on weighted grid using different spatial layers (such as GHSL, LandScan, Nighttime lights) that describes the proportion occupied area within the administrative bounds. To improve the sampling process an additional condition can be provided which ensures that assets are not sampled within a certain radius of preoccupied asset location. Second limitation is the use of simplified correlation between the structural attributes. In the current study only correlation between the building material and building storey is considered in such a way that for masonry type of structure, the maximum number of storey that can be sampled is 4. Such correlation was implemented to prevent generating unrealistic asset. This simplification may be replaced with more sophisticated correlation between the structural attributes (Crowley et al. 2020).

Financial framework for pricing risk-based catastrophe bonds
The financial framework for pricing risk-based cat bonds comprises of three main models: (i) interest rate model; (ii) aggregate loss model; (iii) payoff function. The interest rate model defines the rate of interest to be paid to the bond holder at a specific time of maturity. The Cox-Ingersoll-Ross (CIR) model (Cox et al. 1985) is normally used to price cat bonds thanks to its ability to provide non-negative interest rates. This is based on the assumption of spot interest rate following a mean-reverting square-root process, mathematically expressed by: where k, and are the parameters representing mean-reverting force, long-term mean and the volatility, respectively. W(t) refers to the standard Weiner process for t ∈[0,T]. The condition 2k > 2 ensures that the interest rate remains non-negative. By using the approach proposed by Cox et al. (1985) along with the above-mentioned assumption, the price of a pure-discount T-bond at time t can be computed by where r (t) is the assumed constant to represent market risk.
The aggregate loss model involves computation of aggregate losses, which are defined as the compound distribution of two processes: (i) catastrophe events process, where N(t) denotes the frequency of catastrophic events; (ii) catastrophe severity process, where 1 3 X n defines the severity of losses caused by a catastrophic event. Based on the approach proposed by Ma and Ma (2013), the three following assumptions are usually made: (1) the occurrence of potential catastrophic events of a specific magnitude is modelled as a Poisson point process N (t) (t ∈ [0, T]) , with intensity parameter M . The occurrence time of potential catastrophic events is denoted as (2) the trigger event for the bond is taken at the time when the aggregate loss ( L t ) exceeds a threshold level (D), mathematically expressed as: = inf{t ∶ L(t) ≥ D}; (3) At each time instant ( t n ), the severity of catastrophic events is modelled as an independent and identically distributed (i.i.d) random variable, X n n=1,... , with cumulative distribution function F(x) = P X n < x .
The aggregate loss is defined as: The payoff function models the payment terms for the catastrophe bond contract for the case of trigger and non-trigger events. The payoff function for the zero-coupon bond is given by where Z refers to the face value of bond that is paid to the bond holders when no triggering events occur for the entire duration of contract, that is until maturity time, T. If a triggering event occurs before the maturity time, then the bond holders receive a fraction of the face value. The price of zero-coupon bond is obtained by combining the interest rate model, aggregate model and payoff functions. Under the assumption of risk-neutralized pricing measure Q, the price of zero-coupon ( V ZC ) catastrophe bond at a given time t is mathematically expressed as: where B CIR refers to the spot interest rate at time t; Z is the face value of cat bond (i.e. value of cat bond at time of maturity); F(D, T) denotes the probability function that aggregate losses are less or equal to the loss threshold D, P X 1 + X 2 + ... + X n ≤ D ; m is the intensity parameter of the Poisson point process N(t) at time t; F n (D) is the n-th convolution of F.

Sensitivity analysis
The aim of the sensitivity analysis is twofold. First, it aims to illustrate the implementation of the proposed stochastic exposure model through a realistic example. Second, it aims to quantify the influence of the exposure spatial resolution on the assessment of average annual losses and cat bond prices. To this end, we prepare four exposure models using the proposed approach where the spatial resolution of the models is varied from coarse to fine (i.e. province, municipality and asset level). The sensitivity analysis is carried out on a study region that consists of ten provinces in southern Italy: Avellino, Benevento, Campobasso, Caserta, Foggia, Isernia, Matera, Napoli, Potenza, Salerno. The study area comprises of 879 municipalities, with approximately 1 million (1,186,947) assets. The distribution of assets in each municipality is presented in Fig. 3. We mainly focus on residential buildings, thereby classifying building activity attribute into two sub-groups, residential and others. Figure 4 shows a workflow summary for pricing risk-based catastrophe bonds.
The following sections summarise the development of the hazard, vulnerability and exposure models for the study region and the pricing of the cat bonds considering four exposure models with different spatial resolutions.

Seismic Hazard Model
The seismic hazard model provides the spatial distribution of ground-motion intensity parameter (e.g. peak ground acceleration, PGA), at different return periods (Cornell 1968;McGuire 2004). In this study, we determine the seismic hazard model from a simulation-based seismic hazard analysis (Musson 2000;Crowley and Bommer 2006) that generates a 20,000-year-long synthetic earthquake catalogue, containing 94,446 seismic Area source model for the current study region adopted from SHARE project (Woessner et al. 2015). Grey region refers to the area that is not considered as part of study region (i.e. no events are generated in this region). Admin boundary of selected ten provinces is highlighted in red events. The catalogue is generated by the Monte Carlo sampling approach described by Crowley and Bommer (2006) using the earthquake source model (see Fig. 5) and recurrence parameters provided by the Seismic Hazard Harmonization in Europe (SHARE) project (Woessner et al. 2015). We use the ground motion prediction equation (GMPE) proposed by Bindi et al. (2011) which was calibrated based on the historical strong motion database available from the Italian Accelerometric Archive (ITACA 2011). The inter-event aleatory uncertainty in the GMPE is taken into account by sampling inter from a standard normal distribution for each event in the stochastic catalogue. On the other hand, we neglect the intra-event uncertainty in this study due to the computational expense of the proposed approach. But it is worth noting that intra-event uncertainty accounts for spatial and cross-spatial correlation in the ground motion residual and neglecting such uncertainty may lead to underestimation in the loss estimates (up to 50%) (Silva 2015(Silva , 2018. The local site effects are taken into account through a coefficient in the GMPE that depends on the ground type, as defined by the Eurocode 8 (CEN 2004). To identify the ground type, we consider the shear wave velocity in the upper 30 m V s,30 , shown in Fig. 6. Figure 7 shows the seismic hazard map for a 475-year return period, that is approximately a 10% probability of exceedance in 50 years. It is worth noting that this map incorporates the local site effects using the V s,30 map (see Fig. 6). Figure 8 plots the distribution of PGA within each of the ten provinces. It can be seen that Avellino, Benevento and Potenza have the highest seismic hazard, with a median value of PGA of 0.26 g, 0.24 g and 0.27 g, respectively. The western parts of Caserta and Salerno provinces are characterised by lower seismic hazard, while the north-eastern region is more seismically active as it is adjacent to the seismically active provinces of Benevento and Potenza. The south-western area of Campobasso province that confines to the Benevento province also exhibits relatively high seismic hazard. In the Matera province, the highest seismic hazard is found in the eastern region, where PGA ranges from 0.25  Mori et al. (2020) to 0.37 g. In Foggia, the highest seismic hazard is found in the central and southern parts of the province, where PGA ranges between 0.23 and 0.28 g. The southern area of Napoli province has PGA <0.18 g, while the central and north-western areas have PGA values up to 0.26 g. The seismic hazard of Isernia province is more homogeneous with a value of approximately 0.22 g across the entire province.

Exposure model
Four exposure models, denoted by E1, E2, E3 and E4 (see Table 1), are considered to investigate the influence of the exposure spatial resolution on average annual losses and catastrophe bond pricing. The proposed model is denoted by E1 and is compiled as   . 9 Graphical illustration of random sampling of asset attributes from their respective distribution built based on census data. Dashed line is the representative random number sampled for respective building attribute described in the previous section (see Sect. 2.2). The information used to build exposure model E1 is summarised in Table 2. The asset attributes are sampled from distributions (see Fig. 9) determined from national census data made available by the Italian National Institute of Statistics (ISTAT 2001). In exposure models E2 and E3, assets are relocated and aggregated at the geometric centroid of the municipality (admin 3) and province (admin 2), respectively. The asset attributes for these models are determined from the same distributions used in model E1. Exposure model E4 is built by extracting actual building footprints from OpenStreetMap (OpenStreetMap 2017), which also contains building activity type. The remaining asset attributes: building material, number of storeys and number of dwelling; are sampled from the same distribution used for model E1.

Vulnerability model
The vulnerability model provides the damage ratio as a function of the hazard intensity measure. In the present work, we use the vulnerability curves (see Fig. 10) developed by Rosti et al. (2020) and calibrated based on Italian building damage data collected by Italian Department of Civil Protection (Dolce et al. 2019). To treat the epistemic uncertainty in vulnerability curves, logic trees are implemented for each combination of building material (i.e. Masonry and reinforced concrete) and building storey (L:1-2 storey; MH:>2storeys) as shown in Fig. 11. It is worth to note that for simplicity equal weights are assigned to each branch in the logic tree.

Loss estimation
The losses are computed from the convolution of the hazard, exposure and vulnerability model (presented in previous Sects. 4.1,4.2 and 4.3). For the i th event in the generated stochastic earthquake catalogue, first a stochastic exposure model, SEM i (ref Sect. 2.2) Fig. 10 Vulnerability curves adopted from Rosti et al. (2020) for masonry and reinforced concrete. A:Masonry-irregular layout-flexible floors-with and without tie rods or tie beams, Masonry-irregular layout-rigid floors-without tie rods or tie beams; B:Masonry-irregular layout-rigid floors-with tie beams, Masonry-regular layout-flexible floors-without tie beams; C1:Masonry-regular layout-flexible floors-with tie beams, Masonry-regular layout-rigid floors-with or without tie beams; C2:reinforced concrete buildings designed for both gravity and seismic (pre-1981) loads; D:reinforce concrete building with seismic design post-1981, L:1-2 storey; MH:>2 storey is generated. Thereafter, the loss of each asset in generated SEM i is computed from the damage ratio, floor area (in sqm.) and replacement cost (per sqm.). The damage ratio for an asset is estimated by mapping the PGA value (i.e. PGA experienced by an asset given the i th event) on to their respective vulnerability curves, while the floor area and replacement cost are obtained from the stochastic exposure model. The losses for all assets are then summed together and stored next to the respective event. An event loss table is generated by repeating the above process for all the events in the stochastic earthquake catalogue. It is worth mentioning that for each event in the stochastic earthquake catalogue, a stochastic exposure model is generated which will leverage the Monte Carlo simulation to account for the uncertainty in the location and attributes of the assets in the exposure model. The event loss table is then used to compute the annual probability of exceedance (see Fig. 12) and cumulative distribution function P[L < l] (see Fig. 13). The former is used to calculate average annual losses, defined as the ratio of the sum of losses from each year to the total number of years in the event loss table. The latter is obtained by fitting a lognormal distribution to the simulated loss data from the event loss table. The estimated parameters of the lognormal distribution (i.e. mean and standard deviation) are presented in appendix B. Further, these cumulative distribution functions are used as an input to the aggregate loss model for pricing risk-based catastrophe bonds (see Sect. 4.6).

Influence of exposure spatial resolution on average annual losses
The influence of exposure spatial resolution on average annual loss (AAL) is investigated by considering the model E1 as benchmark case to compare the relative difference in AAL for the model E2, E3 and E4 (see Fig. 14a). Figure 14b shows the relative difference in the footprints extracted from OSM with those extracted from census data. For model E2, AAL estimates are comparable with those computed by the proposed model E1, with a difference ( < 1% ) for most of the provinces, except for Matera, where 8% lower AAL estimates are computed. The E3 model shows higher AAL difference (up to 20% ), possibly because of the uncertainty associated with asset location and low hazard resolution. More specifically such differences are mainly due to the difference in the seismic hazard used in each of the model. This implies that if seismic hazard for E2 and E3 is higher than the E1, then it will lead to higher level of losses as compared to E1 model, vice versa. For example, the seismic hazard at the Fig. 11 Logic trees for each combination of building material and building storey. A:Masonry-irregular layout-flexible floors-with and without tie rods or tie beams, Masonry-irregular layout-rigid floorswithout tie rods or tie beams; B:Masonry-irregular layout-rigid floors-with tie beams, Masonry-regular layout-flexible floors-without tie beams; C1:Masonry-regular layout-flexible floors-with tie beams, Masonry-regular layout-rigid floors-with or without tie beams; C2:reinforced concrete buildings designed for both gravity and seismic (pre-1981) loads; D:reinforce concrete building with seismic design post-1981, L:1-2 storey; MH:>2 storey 1 3 centroid of Matera municipality (E2) and province (E3) is 0.2 g (see Fig. 8) which is lower than median of PGA value (0.23 g) within the entire province; therefore, the AAL values for E2 and E3 are lower than E1 model. The greatest difference in AAL estimates is observed when using the model E4, where the missing footprints results in lower AAL Fig. 12 Comparison of annual probability of exceedance computed from the event loss table for four exposure models: E1 (black solid), E2 (pink dashed), E3 (light blue dashed dot) and E4 (brown dotted). The horizontal red dashed line is at 475-year return period (i.e. 10% probability of exceedance in 50 years). Each point on the curve corresponds to a specific loss value that can be exceeded in given return period values. It is worth noting that more than 50% of assets are missing for the provinces of Avellino, Benevento, Caserta, Isernia and Salerno. On the other hand, for the provinces of Foggia, Matera and Napoli there are more assets than that determined from national census data, suggesting that the OSM data include recently built buildings.
From the above-mentioned observations, we can conclude that E4 behaves poorly for the provinces with insufficient footprint data. Although no significant difference in the annual probability of exceedance is observed between E1 and E2, the difference in AAL 1 3 between these two models cannot be neglected. The model E3 provides higher level of difference in the AAL estimates, possibly because the model has a low exposure and hazard spatial resolutions. Figure 14c presents the distribution of average annual losses for exposure model E1.

Pricing risk-based catastrophe bonds
The framework used for pricing risk-based cat bonds is set out in Sect. 3. We use the 3-month maturity US monthly treasure bills for the period 1994-2013 from Shao et al. (2016) to define the CIR interest rate model. The interest rate model parameters used in Fig. 14 a Relative difference in average annual losses computed for exposure models E2, E3 and E4. using proposed model E1 as benchmark; b relative difference between number of building footprints extracted from OSM and asset locations inferred from national census at province level; c average annual losses for exposure model E1 computed using simulated loss for each municipality. Municipality boundaries within province are denoted by black solid line equations 2-5 are as follows: mean-reverting force k = 0.0984 , long-term mean = 2.04% , volatility = 4.77% , market risk r = −0.01 , initial interest rate r 0 = 2.04%.
For the catastrophe events process, a single-intensity parameter m = 0.3 is used. This is computed as the mean value of the intensity parameters of all seismic source zones in the study region. The cumulative distribution function of losses computed in Sect. 4.4 is used Fig. 15 Contour plot of zero-coupon cat bond price for different loss threshold levels D ∈ [€0.001bn €10bn] and maturity time T ∈ [0.25 3] years to define the mean and standard deviation of the loss function in the catastrophe severity process.
The zero-coupon cat bond is priced with a face value of €1 at time t = 0 . Different times to maturity T and loss threshold levels D are considered: T=[0.25 3] years and D = [€0.001bn €10bn] (see Fig. 15). From Fig. 15, it is observed that for a certain level of time to maturity T, bond price V ZC is directly proportional to the threshold value D, while an indirect relation between bond price and time to maturity is observed for a specific loss threshold value. Moreover, provinces with low probability of exceeding large losses yield higher bond prices. Figure 16 and 17 present the relative difference in bond prices Fig. 16 Relative difference in prices of zero-coupon cat bonds for the following provinces: Avellino, Benevento, Campobasso, Caserta and Foggia. The relative difference is calculated considering the proposed model E1 as benchmark for exposure models E2, E3 and E4 compared to bond price from model E1, our benchmark. No significant difference is observed for prices obtained from model E2, while E4 model yields the largest difference in the bond price. In model E3, higher differences (up to 4.5% ) in bond prices are observed for the larger provinces of Potenza and Salerno, while other smaller provinces display differences ranging between 1 and 3 % . Arguably, these higher differences are observed for larger provinces because the relocation of asset at geometric centroid of province introduces artificial correlation in the asset location and hazard levels. For instance, for the small province of Isernia the difference in bond prices is < 1% . On the other hand, for the larger provinces of Potenza and Salerno, the difference goes up to 6 % . For model E4, the difference in cat bond price depends on missing footprint data, with greater differences observed when more footprint data are lacking. Foggia, Matera and Napoli are the provinces with nearly complete footprint Fig. 17 Relative difference in the price of zero-coupon cat bonds for following provinces: Isernia, Matera, Napoli, Potenza and Salerno. The relative difference is calculated considering the proposed model E1 as benchmark data; thus they show negligible difference < 0.2% . On the other hand, provinces missing significant footprint data, such as Avellino, Benevento, Caserta and Salerno, show a difference ranging between 4 and 9.5%.

Conclusion
The assessment of economic losses caused by earthquakes is paramount for the implementation of effective disaster risk management strategies that use financial instruments, such as insurance, or cat bonds, to improve the financial resilience of communities impacted by strong earthquakes. This study makes an attempt to address the uncertainty in the exposure model by proposing a stochastic model that accounts for the uncertainty in the spatial distributions of assets and their attributes. Furthermore, the paper presents results from a sensitivity analysis carried out to investigate the influence of the spatial resolution of four exposure models on average annual losses (AAL) and cat bond prices. More specifically, in the proposed model (E1), which is used as benchmark, the asset location is generated randomly within each municipality. In model E2, all assets are relocated and aggregated at the centroid of their municipality. In model E3, all assets are relocated and aggregated at the centroid of their province. In model E4, the footprint, and thus the location of each asset, is extracted from publicly available building footprints database. It is worth mentioning that the footprint database can be incomplete, and is noted that a significant number of assets is missing for the provinces of Avellino, Benevento, Caserta, Isernia and Salerno. We conclude that the spatial resolution of the exposure model influences the average annual losses especially for larger urban areas with thousands of assets. More specifically, compared to model E1, model E4 tends to show higher difference of up to 90% , possibly due to the missing building footprints. The AAL difference for model E2 and E3 is 8 % and 20% . On the other hand, cat bond prices computed based on models E1 and E2 are comparable, whereas for models E3 and E4, larger differences are observed (up to 4.5% and 9.5% , respectively). Additionally, we observe that difference in AAL is higher than that observed in cat bond pricing. There are three possible reasons for this, first, the price of the cat bond depends on the loss estimate, interest rate model and pay off function. The second reason is due to the use of broad range of simulated loss data in the cat bond pricing equation where each level of loss is associated with an annual probability of exceedance. On the other hand, AAL is estimated as the average of all these loss values, which may lead to aggregation in the results, thereby exhibiting higher differences compared to cat bond prices. As the current approach has implemented sampling of the asset location from standard uniform distribution, future studies can investigate the use high-resolution weighted grids. These weighted grids can be computed by combining several spatial layers (such as GHSL, LandScan, Nighttime lights), justifying the concentration of population or built areas in the given region. Further improvements can be made by adopting a sophisticated correlation explaining the relationship between asset attributes. Although the proposed stochastic exposure model is developed for seismic hazards, it can be used for other perils where the spatial distribution of assets and their attributes are important for the assessment of the economic losses, specifically from disaster risk management perspective. Table 3 presents recurrence parameters for the source zone used in Sect. 4.1.