3.1.1 Collecting building inventory and calculating exposure

Residential exposure is calculated using building counts by type from the 2016 Canadian Census (Statistics Canada, Government of Canada, 2016) and average square footage estimates from the Canadian Housing Statistics Program (Statistics Canada, Government of Canada, 2020a). For each building type, exposure is computed as:

\begin{equation*}\text{Exposure}=\text{Number of buildings} \times \text{Average area} \times \text{Replacement cost per square foot}.\end{equation*}

Replacement costs are taken from HazCan (Ulmi et al., Reference Ulmi, Wagner, Wojtarowicz, Bancroft, Hastings, Chow, Rivard, Prieto, Journeay, Struik and Nastev2014) and inflated to June 2021 using the Building Construction Price Index (BCPI) (Statistics Canada, Government of Canada, 2020b). To account for regional and estimation uncertainty, we introduce random noise of $\pm 10\%$ to the replacement costs. Content exposure for residential buildings is assumed to be 50% of the structural replacement value (Ulmi et al., Reference Ulmi, Wagner, Wojtarowicz, Bancroft, Hastings, Chow, Rivard, Prieto, Journeay, Struik and Nastev2014; FEMA, 2013). Exposures are then disaggregated by construction material (e.g., wood, steel, concrete) using the general building scheme mapping in HAZUS (Table 5.1 in FEMA (2013)). Exposure is therefore computed separately by building class and construction type, and financial losses are generated using type-specific damage functions applied to the spatially varying shaking intensity. While individual policies are not modeled explicitly, portfolio heterogeneity is captured at the CSD level through building type, construction material, regional variation in floor area, and stochastic variation in replacement costs.

Non-residential exposure is estimated indirectly, as no comprehensive database currently exists. We use building permit data (Statistics Canada, Government of Canada, 2020c) and apply annual ratios of non-residential to residential permits, by type (e.g., commercial, institutional), to infer non-residential square footage. Exposure values for these buildings are then calculated similarly, and construction materials and content assumptions follow the HAZUS framework.

Figure 2 summarizes total calculated exposure per province, with values suppressed to respect the confidentiality of CatIQ’s data, Canada’s Loss and Exposure Indices Provider. See Section 4.1 for details. Western Canada includes British Columbia (BC), Alberta (AB), Saskatchewan (SK), Manitoba (MB), Northwest Territories (NT) and Yukon (YT), and Eastern Canada includes Newfoundland and Labrador (NL), Nova Scotia (NS), Prince Edward Island (PE), New Brunswick (NB), Québec (QC), Ontario (ON), and Nunavut (NU).

Figure 2

Total exposure per province (left) and grouped for Eastern and Western Canada (right). Vertical dashed line separates the two regions.

Figure 3 provides a spatial map of total exposure, building and contents, per square kilometer across CSDs.

Figure 3

Total residential and non-residential exposure (building + contents) per km $^2$ at the CSD level.

Appendix B provides further details on the estimation of exposure.

3.1.2 Simulating earthquake occurrences in space and time

This section briefly summarizes the STPP and residual diagnostics used as supporting tools within the actuarial earthquake loss and capital modeling framework developed in this paper.

A STPP is a stochastic process that models the occurrence of events in both space and time. Applications include the spread of diseases and pandemics, and natural disasters such as earthquakes, tsunamis, and volcanic eruptions. Point processes are studied thoroughly in time; see Cox and Isham (Reference Cox and Isham1980); Daley and Vere-Jones (Reference Vere-Jones2003) and in space; see Cressie (Reference Cressie2015); Moller and Waagepetersen (Reference Moller and Waagepetersen2003); Diggle (Reference Diggle2013). In seismology, STPPs have been widely used to model earthquake patterns, particularly aftershock sequences; see Ogata (Reference Ogata1998); Zhuang et al. (Reference Zhuang, Ogata and Vere-Jones2002).

An STPP is defined on a subset $S\subseteq \mathbb{R}^2\times \mathbb{R}^{+}$ , where each realization consists of space-time points $\lbrace (\boldsymbol{x}_i,t_i) \, : \, i =1, \ldots , n \rbrace$ , where $(\boldsymbol{x}_i,t_i) \in A \times T$ for some known spatial region $A$ and temporal period $T$ , with $\boldsymbol{x}_i\in A\subset \mathbb{R}^{2}$ representing spatial coordinates (longitude and latitude), and $t_i\in T \subset \mathbb{R}^{+}$ the time of occurrence. The first-order properties are described by its spatio-temporal intensity function:

\begin{equation*} \lambda (\boldsymbol{x},t)=\lim _{\lvert d\boldsymbol{x}\rvert ,\lvert dt\rvert \to 0}\frac {\mathbb{E}\left [ N(d\boldsymbol{x},dt) \right ]}{\lvert d\boldsymbol{x}\rvert \lvert dt\rvert }, \end{equation*}

where $d\boldsymbol{x}$ represents a small spatial region around the location $x$ such that $\lvert d\boldsymbol{x}\rvert$ is its area, $\lvert dt\rvert$ represents a small time interval containing the time point $t$ such that $\lvert dt\rvert$ is its length, and $N(d\boldsymbol{x},dt)$ represents the number of events in $d\boldsymbol{x}\times dt$ . Thus, $\lambda (\boldsymbol{x},t)$ represents the expected number of events per unit area per unit of time. If $\lambda (\boldsymbol{x},t)=\lambda$ is constant, the process is said to be homogeneous. Marginal intensities can be defined by integrating over time or space:

\begin{equation*} \lambda _A(\boldsymbol{x})=\int _T\lambda (\boldsymbol{x},t) dt, \quad \text{and} \quad \lambda _X(t)=\int _A\lambda (\boldsymbol{x},t) d\boldsymbol{x}, \end{equation*}

respectively.

Model fit for STPPs can be assessed using the N-test and L-test introduced by Schorlemmer et al. (Reference Schorlemmer, Gerstenberger, Wiemer, Jackson and Rhoades2007). The N-test compares the observed number of events to simulations from the fitted model, while the L-test evaluates whether the observed log-likelihood exceeds most simulated log-likelihoods. These tests evaluate overall fit but provide limited insight into spatial or temporal misfit.

To provide more granular diagnostics, Baddeley et al. (Reference Baddeley, Turner, Møller and Hazelton2005); Zhuang (Reference Zhuang2006) introduced residual analysis methods based on discretizing the domain into a regular grid of pixels. For a fitted intensity $\hat \lambda (\boldsymbol{x},t)$ , the raw residual for a pixel $B_i \in A\times T$ is

\begin{equation*} r^R(B_i) = N(B_i) - \int _{B_i}\hat \lambda (\boldsymbol{x},t) dt d\boldsymbol{x}, \end{equation*}

where $N(B_i)$ is the observed number of events in pixel $B_i$ , which is a measurable set such that $B_i\subset A\times T$ , for $i = 1,\ldots ,m$ , where $m$ is the total number of pre-determined pixels. This compares the total number of observed and estimated points within evenly spaced pixels over a regular grid. Raw residuals can be scaled to improve interpretability. For instance, Pearson residuals rescale the raw residuals to have approximate unit variance, in analogy with Pearson residuals in linear models, i.e.,

\begin{equation*} r^P(B_i) = \sum _{(t_j,\boldsymbol{x}_j)\in B_i} \dfrac {1}{\sqrt {\hat \lambda (\boldsymbol{x}_j,t_j)}} - \int _{B_i}\sqrt {\hat \lambda (\boldsymbol{x},t)}dt d\boldsymbol{x}, \end{equation*}

for all $\hat \lambda (\boldsymbol{x}_i,t_i)\gt 0$ . See Baddeley et al. (Reference Baddeley, Turner, Møller and Hazelton2005) for other methods of scaling the gridded residuals and Baddeley et al. (Reference Baddeley, Møller and Pakes2008) for their properties. Analogous to deviances for regression models, Wong and Schoenberg (Reference Wong and Schoenberg2009) propose using deviance residuals to compare two fitted models $\hat \lambda _1$ and $\hat \lambda _2$ via their log-likelihood contributions:

\begin{align*} r^D(B_i) &= \left (\sum _{(t_j,\boldsymbol{x}_j)\in B_i} \log (\hat \lambda _1(\boldsymbol{x}_j,t_j)) - \int _{B_i} \hat \lambda _1(\boldsymbol{x},t) dt d\boldsymbol{x} \right ) \\ &\quad - \left (\sum _{(t_j,\boldsymbol{x}_j)\in B_i} \log (\hat \lambda _2(\boldsymbol{x}_j,t_j)) - \int _{B_i} \hat \lambda _2(\boldsymbol{x},t) dt d\boldsymbol{x} \right ). \end{align*}

Although pixel-based residuals are informative, they become unstable when the expected number of events per pixel is small, which is a common issue in earthquake modeling. Increasing pixel size to stabilize estimates may distort results in highly inhomogeneous regions.

To address this, Bray et al. (Reference Bray, Wong, Barr and Schoenberg2014) propose Voronoi residual analysis, which partitions the spatial domain using Voronoi tessellations. Each Voronoi cell $C_i$ contains exactly one observed event and comprises all points closer to $\boldsymbol{x}_i$ than to any other point in the process, for all $i=1, 2, \ldots , n$ . Figure 4 illustrates the Voronoi tessellations of earthquake locations in Figure 1. For more details on Voronoi tessellations and their properties, see Boots et al. (Reference Boots, Okabe and Sugihara1999). For a fitted intensity $\widehat \lambda (\boldsymbol{x},t)$ , the raw Voronoi residuals are computed as:

\begin{align*} r^R_V(C_i) &= 1 - \int _{C_i}\hat \lambda (\boldsymbol{x},t) dt d\boldsymbol{x}. \end{align*}

Each residual measures the difference between the observed count (one) and the expected count in cell $C_i$ .

Figure 4

Voronoi tessellation of the earthquake locations, shown in Figure 1.

Building on this, we apply Pearson Voronoi residuals, which scale the residuals in analogy to their pixel-based counterparts, defined by:

\begin{align*} r^P_V(C_i) = \sum _{(t_j,\boldsymbol{x}_j)\in C_i} \dfrac {1}{\sqrt {\hat \lambda (\boldsymbol{x}_j,t_j)}} - \int _{C_i}\sqrt {\hat \lambda (\boldsymbol{x},t)}dt d\boldsymbol{x}. \end{align*}

These residuals are well-defined provided $\hat \lambda (\boldsymbol{x}_i,t_i)\gt 0$ for all $\boldsymbol{x}_i$ and $t_i$ , and yield approximately standard-normal values under model adequacy.

To compare two competing models at a local spatial level, we apply the deviance residual concept to Voronoi polygons following the framework of Bray et al. (Reference Bray, Wong, Barr and Schoenberg2014) and Baddeley et al. (Reference Baddeley, Turner, Møller and Hazelton2005). In analogy with linear models, the resulting residuals may be called deviance Voronoi residuals, defined as such:

(2) \begin{align} r^D_V(C_i) &= \left ( \sum _{(t_j,\boldsymbol{x}_j)\in C_i} \log (\hat \lambda _1(\boldsymbol{x}_j,t_j)) - \int _{C_i} \hat \lambda _1(\boldsymbol{x},t) dt d\boldsymbol{x} \right ) \nonumber \\ &\quad - \left ( \sum _{(t_j,\boldsymbol{x}_j)\in C_i} \log (\hat \lambda _2(\boldsymbol{x}_j,t_j)) - \int _{C_i} \hat \lambda _2(\boldsymbol{x},t) dt d\boldsymbol{x} \right ). \end{align}

Polygons with positive (negative) values indicate regions where $\hat \lambda _1$ ( $\hat \lambda _2$ ) provides a locally improved fit. Summing these residuals across all polygons yields a log-likelihood ratio score that provides a global diagnostic comparison between competing models.

Now, we apply these residual diagnostics to visually and spatially assess the adequacy of candidate spatio-temporal models as a validation step within the broader earthquake loss and capital simulation framework.

To simulate future events, we construct a spatio-temporal point pattern dataset based on the 137 significant earthquakes in Canada between 1900 and 2017 (Lamontagne et al., Reference Lamontagne, Halchuk, Cassidy and Rogers2018). Events prior to 1900 are excluded due to historical incompleteness that could distort the temporal component of the model. We assume space is continuous and time is discrete (annual), and estimate the spatio-temporal intensity function $\hat {\lambda }(\boldsymbol{x}, t)$ using kernel-based methods (Diggle, Reference Diggle1985). This choice is motivated by the well-established robustness and interpretability of kernel smoothing in spatial and STPP modeling. We also follow a practical approach and assume that the STPP is separable such that $\lambda (\boldsymbol{x}, t)$ can be factorized following

(3) \begin{equation} \lambda (\boldsymbol{x}, t)=\lambda _{\boldsymbol{x}}(\boldsymbol{x})\lambda _{T}(t), \end{equation}

for all $(\boldsymbol{x},t) \in A\times T.$ This assumption is justified as our interest lies in rare, high-impact events (e.g., the $\mathrm{PML}_{1/500}$ ), where the exact timing is less critical than the spatial distribution of risk exposure. The spatial component $\lambda _{\boldsymbol{x}}(\boldsymbol{x})$ is the primary focus for Voronoi residual analysis. Model fitting, simulation, and residual computations are conducted in R using the splancs, spatstat, deldir, and stpp packages (Rowlingson and Diggle, Reference Rowlingson and Diggle2022; Baddeley and Turner, Reference Baddeley and Turner2005; Turner, Reference Turner2021; Gabriel et al., Reference Gabriel, Diggle, Rowlingson and Rodriguez-Cortes2022; R Core Team, 2022), along with custom functions developed for Voronoi-based residual diagnostics.

To estimate the temporal component of the intensity function, $\lambda _{T}(t)$ , we apply a Gaussian kernel density estimator with bandwidth $h_T$ selected via Silverman’s rule-of-thumb (Silverman, Reference Silverman1986), i.e.,

\begin{equation*}h_T=0.9\alpha n^{-1/5},\end{equation*}

where $\alpha =\min \lbrace \text{sample standard deviation}, \text{sample interquartile range}/1.34\rbrace$ . For the spatial intensity function $\lambda _{\boldsymbol{x}}(\boldsymbol{x})$ , we use the quartic (biweight) kernel (Diggle, Reference Diggle2013). Bandwidth selection plays a critical role in balancing bias and variance: a smaller bandwidth yields lower bias but higher variance, while a larger bandwidth provides smoother estimates at the cost of increased bias. There is no universally optimal bandwidth, especially for inhomogeneous processes such as earthquakes, which motivates the use of standard data-driven bandwidth selection criteria rather than problem-specific custom estimators. We fit two spatial intensity models using the same kernel but different bandwidths.

• • Model $H_1$ : bandwidth $h_1$ minimizes the estimated mean-square error of $\hat \lambda _{\boldsymbol{x}}(\boldsymbol{x})$ .

• • Model $H_2$ : bandwidth $h_2$ maximizes the leave-one-out likelihood cross-validation score:

  \begin{equation*} LCV(h_2)=\sum _l\log (\lambda _{-l}(\boldsymbol{x}_l))-\int _A \lambda (\boldsymbol{u}) \text d \boldsymbol{u}, \end{equation*}
  where $\lambda _{-l}(\boldsymbol{x}_l)$ is the intensity estimate at $\boldsymbol{x}_l$ excluding that point, and $\lambda (\boldsymbol{u})$ is the smoothed estimate across the domain $A$ . Baddeley et al. (Reference Baddeley, Rubak and Turner2015) discuss the use of likelihood cross-validation in settings where the point pattern exhibits strong clustering, as is typical for earthquake occurrences. Minimizing the mean-square error is used as the primary criterion for global bandwidth selection, while the leave-one-out likelihood-based bandwidth provides an alternative, data-driven smoothing level. Voronoi residuals are subsequently used for spatial diagnostic validation and local comparison of the competing models.

In addition, we fit a homogeneous Poisson process, denoted Model $P$ , with constant intensity $\lambda (\boldsymbol{x})=\lambda$ for all $\boldsymbol{x} \in A$ , to serve as a benchmark for model comparison.

Figure 5

Canada: Raw Voronoi residuals for models $P$ (a), $H_1$ (b), and $H_2$ (c).

Figure 5 displays the raw Voronoi residuals for Models $P$ , $H_1$ and $H_2$ , plotted on the same color scale. An ideal model produces residuals close to zero. Since raw Voronoi residuals are bounded above by 1 (as each cell contains one observed event), large negative values indicate overestimation by the model. Model $P$ clearly underestimates risk in high-seismicity regions of Eastern and Western Canada while overestimating in low-risk zones. This is expected given the homogeneity assumption. Both $H_1$ and $H_2$ perform considerably better, with residuals close to zero across most regions. These visual patterns are confirmed quantitatively by the deviance Voronoi residuals discussed below. Figures 14 and 15 in Appendix C support this conclusion, as do the Pearson Voronoi residual plots in Figures 16–18 therein.

Figure 6

Deviance Voronoi residuals comparing models $H_1$ vs $H_2$ : Canada (a), Western Canada (b), Eastern Canada (c). Positive values indicate regions where model $H_1$ has higher local log-likelihood than $H_2$ .

To further diagnose spatial differences in model fit, we compute deviance Voronoi residuals (see Eq. (2)), plotting the difference in log-likelihood contributions of models $H_1$ and $H_2$ for each Voronoi cell in Figure 6. Model $H_1$ exhibits consistently better local fit than $H_2$ in all regions, particularly in low-risk areas such as Northern and Central Canada, where model $H_2$ tends to overpredict seismicity. In high-risk zones, including the Queen Charlotte fault (west) and the St. Lawrence paleo-rift (east), $H_1$ also performs slightly better. The largest deviance residual occurs in the Queen Charlotte fault zone (Figure 6b). The total log-likelihood ratio score between models $H_1$ and $H_2$ is 210.67, indicating that $H_1$ provides a systematically better local fit than $H_2$ . When compared to model $P$ , the score for $H_1$ is 586.56, while for $H_2$ it is 375.89, indicating that both non-homogeneous models provide substantially improved fit over the homogeneous baseline. Taken together, the mean-square error criterion, likelihood cross-validation, and Voronoi-based residual diagnostics lead us to retain model $H_1$ as the preferred spatial intensity model for the subsequent earthquake loss and capital simulations.

The final fitted spatio-temporal intensity function is obtained by combining the estimated spatial and temporal components:

\begin{equation*}\hat \lambda (\boldsymbol{x}, t)=\hat \lambda _{\boldsymbol{x}}(\boldsymbol{x}) \cdot \hat \lambda _T(t).\end{equation*}

Earthquake simulations are then conducted using the fitted model via the stpp package in R , allowing the generation of a series of synthetic significant earthquake events over multiple years.

3.1.3 Estimating ground shaking intensity

Seismic hazard in Canada is regularly assessed by the GSC, and updates are reflected in revisions to the National Building Code. Ground shaking intensity, a key input in seismic risk modeling, is typically expressed through attenuation relationships of ground motion indices such as peak ground acceleration (PGA), peak ground velocity, or spectral acceleration. This component represents the physics-informed layer of our otherwise data-driven simulation framework, allowing physical ground-motion processes to directly inform the actuarial loss and capital calculations.

In the proposed framework, earthquake magnitudes are not simulated from a standalone parametric magnitude distribution. Instead, magnitude enters implicitly through the physics-informed ground-shaking component. For each simulated event location, PGA is sampled from a Generalized Pareto Distribution (GPD) fitted to the discrete GSC seismic hazard exceedance points at the nearest grid location. PGA is then converted to MMI using the empirical relationship of Wald et al. (Reference Wald, Quitoriano, Heaton and Kanamori1999). The spatial decay of shaking with distance from the epicenter is modeled using the attenuation relationships of Bakun et al. (Reference Bakun, Johnston and Hopper2003) and Bakun and Wentworth (Reference Bakun and Wentworth1997), which provide an explicit functional link between intensity, magnitude, and source distance.

We use the GSC’s seismic hazard values for firm soil sites from the 2015 National Building Code of Canada (NBCC), covering a grid over Canada and adjacent areas (Halchuk et al., Reference Halchuk, Adams and Allen2015). For each grid point, the mean PGA is available at eight annual exceedance probabilities: 0.02, 0.01375, 0.0100, 0.00445, 0.0021, 0.0010, 0.0005, and 0.000404. A Generalized Pareto Distribution (GPD) is fitted at each grid point to interpolate PGA values across the various return periods. The quantile function for the GPD is provided in Eq. (9) in Appendix D.

For each simulated earthquake epicenter, we identify the nearest PGA grid point and sample a ground motion value using its fitted GPD.

The corresponding MMI level is then estimated using the empirical relationship:

\begin{equation*}\text{MMI} = 3.66 \log _{10}(\text{PGA}) -1.66,\end{equation*}

with standard error 1.08, for PGA in cm/s $^2$ (Wald et al., Reference Wald, Quitoriano, Heaton and Kanamori1999).

To model the spatial distribution of shaking, we use attenuation equations developed by Bakun et al. (Reference Bakun, Johnston and Hopper2003) and Bakun and Wentworth (Reference Bakun and Wentworth1997) to estimate moment magnitude as a function of MMI and epicentral distance d (in km):

• • Eastern Canada (Bakun et al., Reference Bakun, Johnston and Hopper2003):

  (4) \begin{equation} \text{M} = \frac {1}{1.68}\left (\text{MMI} - 1.41 + 0.00345\text{d} + 2.08\log _{10}(\text{d})\right ), \end{equation}

• • Western Canada (Bakun and Wentworth, Reference Bakun and Wentworth1997):

  (5) \begin{equation} \text{M} = \frac {1}{1.09}\left (\text{MMI} - 5.07 + 3.69\log _{10}(\text{d})\right ), \end{equation}

To ensure consistency with the definition of a significant earthquake (Lamontagne et al., Reference Lamontagne, Halchuk, Cassidy and Rogers2018), only simulated events whose inferred magnitudes, obtained after converting PGA to MMI and applying the attenuation relationships above, exceed $M\gt 6$ are retained for the loss analysis. Since a single earthquake can generate multiple intensity values at varying distances, we construct isoseismal maps, which are contours of equal MMI, by inverting the above relationships to determine the distance d corresponding to each MMI level. Figure 7 shows a sample isoseismal map produced for a simulated earthquake. The concentric shaded regions correspond to different MMI thresholds, illustrating how intensity attenuates with distance from the epicenter of the earthquake. We account for parameter uncertainty by introducing random noise based on the standard errors reported in Bakun and Wentworth (Reference Bakun and Wentworth1997); Bakun et al. (Reference Bakun, Johnston and Hopper2003); Wald et al. (Reference Wald, Quitoriano, Heaton and Kanamori1999). For each MMI level, we identify the affected CSDs and calculate the proportion of each CSD’s area falling within the corresponding isoseismal circle. Algorithm 1 outlines the procedure for simulating an earthquake and generating the corresponding isoseismal maps.

Figure 7

Isoseismal map for a simulated earthquake, showing the spatial extent of MMI VI–X.

3.1.4 Estimating earthquake-induced damage rates

Earthquakes may result in both financial and non-financial losses. This study focuses on direct financial losses due to physical damage to buildings and their contents. Losses from secondary hazards (e.g., landslides, tsunamis, fire) and non-financial losses (e.g., fatalities, business interruption) are outside the scope of this analysis.

Structural and non-structural damage are modeled separately:

• • Structural damage refers to the building’s load-bearing components (e.g., walls, roofs).

• • Non-structural damage is divided into:

  1. – Drift-sensitive components (e.g., partitions, cladding),

  2. – Acceleration-sensitive components (e.g., mechanical systems, elevators),

  3. – Building contents (e.g., furniture, electronics).

To quantify damage, we rely on MMI-based damage probability matrices (DPMs), which assign to each building type the probability of being in one of seven defined damage states (e.g., no damage, slight, moderate, destroyed, $\ldots$ ) at each MMI level. DPMs were developed by Rojahn et al. (Reference Rojahn and Sharpe1985) and extended to BC by Thibert (Reference Thibert2008) under assumptions of firm ground, regular geometry, and pre-1990 building codes.

Table 1 shows an illustrative example of a DPM, presenting probabilities by damage state for wood-frame residential buildings at MMI levels VI to XII. Each damage state is associated with a damage factor range, which is the proportion of building value expected to be lost.

Table 1.

DPM for structural damage in wood light frame residential building (Thibert, Reference Thibert2008)

^*** represents very small probability values, almost 0.

To quantify damage as a proportion of exposure, we compute the mean damage factor (MDF) for each CSD, building type, and damage component. The MDF is the expected damage percentage at a given MMI level and is computed as the weighted sum of damage factors based on their associated probabilities in the DPM. To introduce simulation variability, damage factors are randomly sampled from their corresponding ranges. Separate MDFs are computed for each damage type (structural (S), acceleration-sensitive (AS) non-structural, drift-sensitive (DS) non-structural, and building contents (BldgC)). These are computed for each CSD-MMI pair and subsequently used to estimate insured losses in Section 3.1(v).

3.1.5 Estimating seismic losses and insurance claim payments

For each simulated event, component losses in CSD $j$ and building class $k$ are computed by applying mean damage factors (MDF) to exposures allocated among structural (S), drift-sensitive (DS), acceleration-sensitive (AS), and contents (BldgC) components. Following Onur et al. (Reference Onur, Ventura and Finn2005), exposures are split as 25% to S, 37.5% to DS, and 37.5% to AS; contents exposure is treated separately. Denoting

\begin{equation*} E_{j,k}=\text{building exposure}_{j,k}, \quad C_{j,k}=\text{contents exposure}_{j,k},\end{equation*}

the component losses at MMI level $m$ are

\begin{align*} L_{j,k,S}^{(m)} &= \mathrm{MDF}_{S,k}^{(m)} \times 0.25\,E_{j,k},\\[3pt] L_{j,k,DS}^{(m)} &= \mathrm{MDF}_{DS,k}^{(m)} \times 0.375\,E_{j,k},\\[3pt] L_{j,k,AS}^{(m)} &= \mathrm{MDF}_{AS,k}^{(m)} \times 0.375\,E_{j,k},\\[3pt] L_{j,k,BldgC}^{(m)}&= \mathrm{MDF}_{BldgC,k}^{(m)} \times C_{j,k}, \end{align*}

and the total loss is

\begin{equation*} \mathrm{Loss}_{j,k}^{(m)} = L_{j,k,S}^{(m)} + L_{j,k,DS}^{(m)} + L_{j,k,AS}^{(m)} + L_{j,k,BldgC}^{(m)}. \end{equation*}

When a CSD spans multiple MMI levels, each $L_{j,k,*}^{(m)}$ is weighted by the area fraction in level $m$ .

Table 2 shows market penetration $\pi _j$ , deductible $d_j$ , and policy limit $u_j$ used in insurance claim payment calculations. For each CSD:

\begin{align*} \mathrm{Claim}_{j,k} = \begin{cases} 0, & \mathrm{Loss}_{j,k}\le d_j,\\[3pt] \pi _j\,(\mathrm{Loss}_{j,k}-d_j), & d_j\lt \mathrm{Loss}_{j,k}\le u_j,\\[3pt] \pi _j\,(u_j-d_j), & \mathrm{Loss}_{j,k}\gt u_j, \end{cases} \end{align*}

where

\begin{equation*} d_j = \delta _j \sum _{k} E_{j,k}, \quad u_j = \lambda _j \sum _{k} E_{j,k}, \end{equation*}

with $\delta _j$ and $\lambda _j$ the deductible and limit percentages at CSD $j$ . Total losses and claims per CSD are

Table 2.

Deductibles, policy limits, and market penetration rates (AIR Worldwide, 2013)

Penetration = fraction of exposed value insured. Deductible and limit = percentage of insured exposure.

\begin{equation*} \mathrm{Loss}_j = \sum _k \mathrm{Loss}_{j,k}, \quad \mathrm{Claim}_j = \sum _k \mathrm{Claim}_{j,k}. \end{equation*}

Appendix E details Algorithms 1 and 2, which execute the simulation of earthquake events and the calculation of losses and claims. Together, Sections 3.1(ii)–3.1(v) form an end-to-end, fully reproducible actuarial pipeline that links spatio-temporal earthquake occurrence, physics-informed ground shaking, exposure by building type, and insurance loss and claim modeling within a single simulation framework.