# Bayesian modeling of insurance claims for hail damage

Ophélia<EMAIL_ADDRESS>[ Anthony C.
<EMAIL_ADDRESS>[ Timo
<EMAIL_ADDRESS>[ Institute of Mathematics, École
Polytechnique Fédérale de Lausanne (EPFL), Switzerland  Institute for
Environmental Decisions, ETH Zurich, Switzerland  Federal Office of
Meteorology and Climatology MeteoSwiss, Zurich, Switzerland
presep= ]e1presep=, ]e2presep=, ]e3

###### Abstract

Despite its importance for insurance, there is almost no literature on
statistical hail damage modeling. Statistical models for hailstorms exist,
though they are generally not open-source, but no study appears to have
developed a stochastic hail impact function. In this paper, we use hail-
related insurance claim data to build a Gaussian line process with extreme
marks to model both the geographical footprint of a hailstorm and the damage
to buildings that hailstones can cause. We build a model for the claim counts
and claim values, and compare it to the use of a benchmark deterministic hail
impact function. Our model proves to be better than the benchmark at capturing
hail spatial patterns and allows for localized and extreme damage, which is
seen in the insurance data. The evaluation of both the claim counts and value
predictions shows that performance is improved compared to the benchmark,
especially for extreme damage. Our model appears to be the first to provide
realistic estimates for hail damage to individual buildings.

Climate extremes,

Extreme value analysis,

Hail damage,

Line process,

Bayesian modeling,

###### keywords:

, , ,

##  1 Introduction

Global warming has already begun to affect the behaviour of insurers
worldwide, both by increasing premiums and by making companies unwilling to
underwrite some risks; a recent example is the May 2023 decision by the US
company State Farm to cease offering house insurance to new clients in
California. Hail is of particular interest to Swiss insurance companies
because of large annual insured losses, averaging several million Swiss francs
(CHF) (Botzen, Bouwer and van den Bergh, 2010), though there is substantial
year-to-year variability. Northern Switzerland experienced significant hail
events in 2021, leading to estimated insured losses of 2 billion CHF and
placing a heavy financial burden on insurers (Müller, 2021). The risk of large
losses due to hailstorms is increasing as a result of the construction of more
new buildings each year, and climate change may increase the frequency of
damaging hailstorms in Europe more broadly (Rädler et al., 2019). Destructive
hailstorms are also an important risk for agriculture, buildings, and vehicles
elsewhere in the world (Changnon, 2008; Warren et al., 2020), but despite
their importance, such storms can be very localised and are hard to model.

The literature on the statistical modeling of the impact of hailstorms on
buildings is very limited. Although stochastic models for hailstorm risk
(Deepen, 2006; Otto, 2009; Punge et al., 2014; Púčik et al., 2017) or
hailstone size (Perera et al., 2018; Liu, Li and Wang, 2021) exist, most open-
source studies on hailstorm impact use deterministic functions to link the
intensity of a hail hazard to its monetary damage. The spatial footprint of
hail events has been discussed in a few recent modeling studies, in which
hailstorms are either represented as ellipses (Otto, 2009) or stretches of
constant width (Deepen, 2006). Punge et al. (2014) use a Poisson distribution
to estimate the frequency of hail in Europe on a $50\times 30$km grid, using a
bimodal normal distribution in each grid cell to estimate the pointwise
probability of hail based on hail reports and observations of overshooting
cloud tops, whose presence for over ten minutes can indicate thunderstorm
severity. Deepen (2006) generates a stochastic catalog of hail events in
Germany by simulating areas of fixed width, random location, and random
length. Hailstorms are also simulated by modeling random hailstones in Liu, Li
and Wang (2021), while Perera et al. (2018) propose a hailstone size
distribution to aid in impact estimation.

On the damage modeling side, Hohl, Schiesser and Aller (2002) propose a
logistic impact function derived from hail kinetic energy, and Schuster, Blong
and McAneney (2006) further explore the link between this energy and impacts.
Claim data from insurance companies, though not usually publicly accessible,
is a valuable source of information on hail impacts and has been used to
complete the radar signal for hail in several recent studies. The hail damage
model developed in Schmidberger, Daniell and Kunz (2018), for example, derives
hail tracks in Germany from radar and insurance data, and Deepen (2006) uses
simulated hail footprints to model damage to cars through a Poisson
distribution fitted with vehicle insurance data. Brown, Pogorzelski and
Giammanco (2015) use insurance data to explore the link between roof material
and hail impact on buildings in Texas.

These studies all involve randomness from the hail event itself. Indeed,
radar-based proxies are used to derive the probability and/or the expected
intensity of a hailstorm on grids with resolution of several kilometers. Those
proxies are often chosen over direct hail measurements with automatic hail
sensors (Kopp et al., 2023) because they are more spatially consistent.
However, the monetary impacts due to hailstorms appear in very narrow and
localized tracks that are usually poorly represented by models with such
grids.

The goal of the present study is to propose a spatially consistent model for
insurance claims related to hail damage at the building level. This model
differs from previous ones, as the probability and intensity of a hail event
are supposed to be known, and stochasticity comes from the possible spatial
impacts of a hail storm. The model we develop seems to be the first to combine
a random line process and an extreme value model in order to represent hail
damage tracks accurately. We describe the data available to us in Section 2,
introduce a line model with extreme marks in Section 3, and describe the
results we obtain when applying this model in the Swiss canton of Zürich,
henceforth “the canton”, in Section 6.

##  2 Data and initial analysis

###  2.1 Data

Two variables representing hail risk were provided by the meteorological
service MeteoSwiss in the scope of the scClim project, the purpose of which is
to combine knowledge from different fields to create a continuous model chain
from simulating thunderstorms to quantifying the monetary impacts of hail in
Switzerland. Gridded one-kilometer-resolution maps of the probability of hail
(POH) and the maximum expected severe hail size (MESHS), derived from
volumetric radar reflectivity (Nisi et al., 2016), are available for the
canton during the convective season (April–September) for the years 2002–2021.
The MESHS offers more spatial granularity than the POH and thus was preferred
(Figure 1).

The output from a MESHS-based deterministic damage function developed during
the scClim project (Schmid et al., 2023+) and calibrated with insurance data
in the canton is also available for the same period; the results in that paper
use the same function, calibrated over several Swiss cantons. The data made
available to us contain predictions for the numbers of affected assets and the
monetary hail damage in each cell of a 2km square grid. For conciseness below
we shall use the terms “grid cell” or sometimes just “cell” in reference to
this grid. The PAA and damage functions were developed following the
hazard/exposure/vulnerability methodology of the CLIMADA framework described
by Aznar-Siguan and Bresch (2019) and will be simply referred to as “CLIMADA”
below. We use a per-building version of the CLIMADA output, referred to as
“downscaled CLIMADA” below, that will be used as input for our claim value
model. This per-building damage is a naive downscaling of the per-cell CLIMADA
damage, and attributes weights to each building as a function of its insured
value. In practice, this means that every building in the cell is impacted
when CLIMADA predicts a positive per-cell value, artificially inflating the
number of buildings affected.

In addition to hail-related variables, we use wind direction from the state-
of-the-art ERA5 reanalysis from the European Centre for Medium-Range Weather
Forecasting (Hersbach et al., 2020), available from 1979 onwards on a 25km
square grid over Europe.

Figure 1: Hail risk covariates POH (left), MESHS (center), and wind direction
(right) on 28 June 2021 in the canton ($47.15$–$47.70^{\circ}$N,
$8.35$–$8.99^{\circ}$E).

Insurance data for hail-related claim damage to buildings in the canton is
also available from one of the stakeholders of scClim, the Zürich cantonal
insurance company GVZ. These data consist of individual claim values for
buildings for the period 2000–2022, during which there were 244 days with
positive claims somewhere in the canton and a total of 46254 claims. These are
the amount finally paid by the insurance company in Swiss francs (CHF)
following hail damage to a building and not estimated values of monetary
damage. The construction year, volume, and actualized insured value of every
insured building in the canton are also available. We did not explore
potential issues linked with preferential sampling, since the owner of every
building in the canton is legally obliged to take out natural hazards
insurance with GVZ. Consequently, data are very dense, though there are
spatial disparities in exposure owing to variations in population density. As
exposure equals the insured monetary value of buildings, urban areas are much
more exposed than suburban areas or the countryside, but the distribution of
buildings is spatially rather homogeneous in terms of average exposed value
per cubic meter; see Figure 2. Our claim value model should be able to predict
the difference between the true damage, i.e., the claim values reported by
GVZ, and downscaled CLIMADA damage when a positive claim was recorded by GVZ.
We call this target variable the “residual damage”.

Figure 2: Insured monetary value for individual buildings (left) and the
exposed monetary value per $m^{3}$ (right).

###  2.2 Exploratory analysis

####  2.2.1 Hail footprints

Literature exploring hailstorm patterns agrees on an ellipsoidal shape (Otto,
2009; Punge et al., 2014), or a sufficiently wide straight sketch (Deepen,
2006) for modeling the spatial extent of a hail footprint. Two of those
studies explore hail tracks in Germany (Deepen, 2006; Otto, 2009) while the
third treats hailstorm footprints in Central Europe (Punge et al., 2014).
Although modeling single hail events on a large territory with locally bounded
shapes seems reasonable, the canton of Zürich covers a much smaller area than
Germany or Central Europe, and individual hail-related claims in our data
suggest that hail storms progress along a line in space with a very narrow
lateral dispersion; see Figure 3. Exploratory analysis suggested that the line
direction is related to the average wind direction on the day of a hail storm.

Figure 3: Example days with more than 50 recorded claims. Red-colored squares
in the first row correspond to the aggregated claim count per 2km grid cell.
The red line has a slope corresponding to the average wind direction on that
day. The graphs on the second row represent the distance from the centroid of
a 2km grid cell to the red line.

####  2.2.2 Local and extreme damage

The deterministic damage function developed through CLIMADA provides good
estimates for the amount of monetary damage on a grid and adequately
represents spatial patterns over the canton. However, the naively downscaled
CLIMADA compensates mispredicted individual claim damage values with a very
large number of claims (see Figure 4) as damage is distributed over all
exposed buildings within a cell. Indeed, CLIMADA cannot distinguish between a
few claims of high damage and many claims with low damage in a cell.
Furthermore, increases in the frequency or intensity of hail would impact
either the count or the value of hail-related damage, causing the compensation
mechanism described above to fail. There is thus a need for a model to provide
realistic values of the damage per building, which is highly relevant for
insurance. The objective of this study is to provide such a model for hail-
related monetary damage that respects the count/size ratio observed in the
claim dat, by lowering the frequency of positive claims while allowing their
values to be locally extreme.

Figure 4: Number of buildings with positive predicted damage per grid cell per
day using the downscaled CLIMADA impact function: (a) claim value per single
building per day: (b) and total damage recorded on a given day aggregated over
the canton: (c) all with the line $x=y$.

Figure 5 shows strong seasonal variation in claims: most damage occurs between
June and August, with a peak in June, and claims occurring in April or
September look less heavy-tailed than in May–August. Henceforth we only
consider claims in April–September, which represent 99.74% of the total data,
and define the “hail season” to be May–August.

Figure 5: Exploratory analysis of claim values: (a) average total claim value
over the canton per month: (b) Spearman correlation $\rho$ and extremal
correlation $\pi$ for pairwise time series of claim values per cell as a
function of the distance between cells. In (b) the solid lines represent the
average pairwise correlation over equally distant cells, and the shaded areas
show the $90\%$ confidence range.

###  2.3 Spatiotemporal correlation

Claims for damage from a hail event can be made on that day or with a lag of a
few days, so the damage function derived through CLIMADA pre-processes the
original insurance data to cluster claims received during a 4-day window
around a big hail event detected with the POH values. For any reported claim,
POH values in a $\pm 2$-day window are scanned. If a POH higher than 50% of
that on the day of the reported claim is observed, the claim date is changed
to the day with the highest POH (Schmid et al., 2023+). This pre-processing
step largely succeeds in removing short-term autocorrelation in the claim
values, so here we focus on spatial correlation. For this purpose, we
introduce the extremal correlation $\pi_{h}(u)=\mathbb{P}(X_{s+h}\geq u\mid
X_{s}\geq u)$ of the variable $X_{s}$, where $s$ denotes the spatial location
of a cell, $u$ is a high threshold and $h$ a spatial lag. The threshold $u$ is
chosen by applying the threshold selection method described in Varty et al.
(2021) to the log total sum of damage $X_{s}$ (see Supplementary Material). We
also study the pairwise Spearman correlation $\rho$ for the daily sum of claim
values per cell. Figure 5 shows that both $\pi$ and $\rho$ decrease as the
distance between two spatial locations increases.

##  3 Key model elements

In this section we describe the key elements of a model to reproduce the very
localized but large damage seen in the data. We first explain how the long and
narrow hail footprint (as observed in Section 2.2.1) can be modeled using a
Gaussian line process, and then recall the peaks-over-threshold approach from
extreme value theory, which is used to account for large impacts on individual
buildings (see Section 2.2.2). More details of the modeling are given in
Section 4.

###  3.1 Random line process

Claim values are usually represented as a spatiotemporal point process
$s_{t}=(t,x,y)$, with $t\geq 0$ and $(x,y)$ the geographical coordinates. In
the forthcoming discussion, a sequence of distance-conserving transformations
will be applied to map this onto a coordinate system that is better suited for
defining the random line model.

For $t>0$ let $\Theta_{t}\in[-\pi,\pi]$ and $\alpha_{t}$ respectively be the
time-varying random angle and vertical deviation of a line $\mathcal{L}$ from
a chosen origin point $s_{0}=(x_{0},y_{0})$, defined as the set of points

$\mathcal{L}=\left\\{(s,t):s^{\mathrm{\scriptscriptstyle
T}}\begin{pmatrix}-\tan{\Theta_{t}}\\\ 1\end{pmatrix}=\alpha_{t}\right\\},$
(1)

where $s=(x,y)$ denotes geographical coordinates and $t$ is the discrete time
coordinate. At time $t$ the projection of any pair of spatial coordinates
$s=(x,y)$ in the coordinate system in which $\mathcal{L}$ is the horizontal
axis and the origin is $s_{0}$ may be written as

$s_{t}^{\mathcal{L}}=\begin{pmatrix}\cos{\Theta_{t}}&\sin{\Theta_{t}}\\\
-\sin{\Theta_{t}}&\cos{\Theta_{t}}\end{pmatrix}\begin{pmatrix}x\\\
y-\alpha_{t}\end{pmatrix},$ (2)

and the orthogonal projection of $s$ onto the line would thus be the point
$\pi_{t}^{\mathcal{L}}(s)=(\cos{\Theta_{t}x}+\sin{\Theta_{t}}y,0)$ in the new
coordinate system. The Euclidean distance between a point in space and the
line $\mathcal{L}$ can be computed in any coordinate system, and as
$s_{t}^{\mathcal{L}}$ and $\pi_{t}^{\mathcal{L}}(s)$ have the same
$x$-coordinate in the new system, at time $t$ this distance can be expressed
as

$d_{t}\left(s,\mathcal{L}\right)=|y^{\prime}|=|(y-\alpha_{t})\cos{\Theta_{t}}-x\sin{\Theta_{t}}|.$
(3)

To allow the intensity of points at time $t$ to be highest close to the random
line $\mathcal{L}$, we define a spatiotemporal Gaussian field $X^{\mu}(s,t)$
whose mean is

$m_{t}(s)=\dfrac{\sigma_{m}}{1+d_{t}\left(s,\mathcal{L}\right)}-1,$ (4)

where $\sigma_{m}$ is a dispersion parameter that controls the concentration
of points around $\mathcal{L}$. In accordance with the exploratory analysis,
we chose a correlation function $\rho$ such that the correlation between
$X^{\mu}(s_{0},t)$ and $X^{\mu}(s_{1},t)$ decreases when the distance $w$
between $s_{0}$ and $s_{1}$ increases (see Figure 5). A common choice is the
Matérn correlation function

$\rho(w)=\\{2^{\nu-1}\Gamma(\nu)\\}^{-1}(w/l)^{\nu}K_{\nu}(w/l),\quad w>0,$
(5)

where $\nu>0$ is a shape parameter controlling the smoothness of the Gaussian
process, $l>0$ is a scale parameter, $\Gamma(\cdot)$ is the Gamma function,
and $K_{\nu}(\cdot)$ is the modified Bessel function of the second kind. After
some experimentation, we took $\nu=1.5$, which gives fields of similar
smoothness to the data, and estimate the parameter $l$ as part of a
hierarchical Bayesian model. When $\nu=1.5$, Equation (5) simplifies to

$\rho(w)=\left(1+{\sqrt{3}w}/{l}\right)\exp{\left(-{\sqrt{3}w}/{l}\right)},\quad
w>0.$ (6)

###  3.2 Marginal model for extreme claim values

The exploratory analysis suggests using extreme value theory to model the
largest claim values. The generalized Pareto distribution,

$\displaystyle{\rm GPD}^{u}(x)$
$\displaystyle=1-\left(1+\xi\dfrac{x-u}{\sigma_{u}}\right)_{+}^{-1/\xi},\quad
x\in[u,\infty),$ (7)

where $a_{+}=\max(a,0)$ for real numbers $a$, provides a standard model for
the exceedances of a high threshold $u$. The model depends on a shape
parameter $\xi$ that determines the weight of the distribution tails and on a
scale parameter $\sigma_{u}$; both are specified in Section 4.2.2. We select a
constant threshold $u$ by applying the method described by Varty et al. (2021)
to the log of the total sum of claim values over the canton; see the
Supplementary Material.

##  4 Modeling extreme hailstorms

Our model for hail damage uses a discrete zero-inflated count process for the
number of claims and a continuous two-part distribution for hail damage
values. In the following section, the spatiotemporal matrix of observed
covariates will be designated by the letter $M$. If a variable needs to be
specified, it is written $M^{\text{NAME}}$ — for the MESHS, for instance, we
write $M^{\text{MESHS}}$. Building exposure, MESHS, POH, CLIMADA predicted
claim count and downscaled CLIMADA value are respectively designated by Exp,
MESHS, POH, NC, and YC. The count of individual claims in a grid cell on a
given day is denoted by $N$, while the value of an individual claim is
designated by $Y$.

###  4.1 Hail damage count

The claim count is modeled on a 2km square grid. Hail can either strike very
locally and violently or can be spread out more smoothly, as observed in the
data, in which the maximum number of observed claims in a cell on a single day
is 469, and the nonzero minimum is 1. In view of this wide range of values, we
model $N$ as a negative binomial. Positive counts are scarce, so we model them
as realisations of a zero-inflated negative binomial random variable, with
probability mass function

$\text{NB}_{\psi,\mu,\alpha}(x)=\begin{cases}(1-\psi)+\psi\left(\dfrac{\alpha}{\alpha+\mu}\right)^{\alpha},&x=0,\\\
\psi\dfrac{\Gamma(x+\alpha)}{x!\Gamma(\alpha)}\left(\dfrac{\alpha}{\mu+\alpha}\right)^{\alpha}\left(\dfrac{\mu}{\mu+\alpha}\right)^{x},&x=1,2,3,\ldots,\end{cases}$
(8)

where $\psi\in(0,1)$, $\mu>0$ and $\alpha>0$ is a shape parameter. We set
$N\mid\psi,\mu,\alpha\sim\text{NB}_{\psi,{\mu},\alpha}$.

The probability $\psi$ of observing a non-zero claim in grid cell $s$ on day
$t$ is modeled as

$\psi(s,t)={\rm
expit}\left\\{\psi_{0}+\psi_{1}\mathbf{1}_{M^{\text{NC}}_{s,t}>0}+\psi_{2}M^{\text{NC}}_{s,t}m_{t}(s)\right\\},$
(9)

where ${\rm expit}(x)=\\{1+\exp(-x)\\}^{-1}$ and $\psi_{0},\psi_{1},\ldots$
are real parameters. We define the mean $\mu$ of the negative binomial
variable through the equation

$\displaystyle\log\mu(s,t)=\mu_{0}+\sum_{i=1}^{3}\mu_{1i}\left(M_{s,t}^{\text{NC}}\right)^{i}+\mu_{2}M_{s,t}^{\text{NC}}m_{t}(s)+X^{\mu}(s,t)+\epsilon(t),$
(10)

where $X^{\mu}$ is a spatiotemporal Gaussian field whose mean $m_{t}$ and
covariance function $\rho$ are given respectively in (4) and (6), and
$\mu_{0},\ldots$ are real parameters. The Gaussian noise $\epsilon(t)$ has
mean zero and one variance for the months April and September and another
variance for the months May–August.

###  4.2 Hail damage values

The number of positive claims (46,254) is much smaller than the roughly
350,000 cell-date combinations in which hail events might have occurred, so it
is reasonable to model spatial patterns at a coarser resolution than the 2km
grid used for the counts. Spatial Gaussian random fields used to model
unobserved covariates underlying the hail damage values are thus defined over
a grid of resolution roughly 10km in which each cell has at least 100 positive
claims over the years 2000–2015.

Downscaled CLIMADA under-predicts the values of $99.3\%$ of reported claims,
so we model only positive errors, i.e., if the predicted count $N$ is
positive, we only allow a shift upwards from the downscaled CLIMADA value
$M^{\text{YC}}$. The resulting residual hail damage variable
$Z=Y-M^{\text{YC}}$ is modelled using a beta model for non-extreme values and
a generalised Pareto model for extreme values.

We first introduce a binary variable $R$ to model the event that a claim value
exceeds the threshold $u$ ($R=1$) or not ($R=0$), with success probability

$p(s,t)={\rm
expit}\left\\{p_{0}+p_{1}M_{s,t}^{\text{POH}}+p_{2}M_{s,t}^{\text{MESHS}}+p_{3}M_{s,t}^{\text{MESHS
$\cdot$ POH}}+p_{4}M_{s,t}^{\text{Exp}}+\chi(s)+\epsilon_{p}(t)\right\\},$
(11)

where $M_{s,t}^{\text{MESHS $\cdot$
POH}}=M_{s,t}^{\text{MESHS}}M_{s,t}^{\text{POH}}$ and $\chi$ and
$\epsilon_{p}$ normally-distributed random effects respectively per grid cell
and season. The Beta and Pareto models for non-extreme and extreme claim
values are detailed below. In the following sections, $f$ denotes the function
$f:x\mapsto\log(1+x)$.

####  4.2.1 Non-extreme residual damage

Residual damage $Z$ for which $f(Z)\leq u$ is described by letting
${Z}/{f^{-1}(u)}$ have a beta density with mean $\nu$ and variance
${\nu(1-\nu)}/{\kappa+1}$,

$\textrm{Beta}_{\nu,\kappa}(x)=\dfrac{x^{\nu\kappa-1}(1-x)^{(1-\nu)\kappa-1}}{B\left\\{\nu\kappa,(1-\nu)\kappa\right\\}},\quad
x\in(0,1),$ (12)

where $B(\alpha,\beta)={\Gamma(\alpha)\Gamma(\beta)}/{\Gamma(\alpha+\beta)}$.
We set

$\dfrac{Z}{f^{-1}(u)}\quad\mid\quad\left\\{f(Z)\leq
u\right\\},\mu_{B},\sigma_{B}\sim\text{Beta}_{\nu,\kappa},$ (13)

and model the mean of this variable via the expression

$\nu_{t}(s)={\rm
expit}\left\\{\nu_{0}+\nu_{1}M_{s,t}^{\text{POH}}+\nu_{2}M_{s,t}^{\text{MESHS}}+\nu_{3}M_{s}^{\text{Exp}}+X^{\beta}(s)\right\\},$
(14)

where $X^{\beta}$ is a spatial Gaussian process with zero mean and covariance
kernel function

$\rho_{\alpha}(w)=\left(1+\dfrac{w}{4l^{2}}\right)^{-2}.$ (15)

####  4.2.2 Extreme residual damage

Damage arising when $f(Z)>u$ is modeled by letting $f(Z)-u$ be a generalized
Pareto variable,

$f(Z)-u\,\mid\,\\{f(Z)>u\\},\sigma_{u},\xi\sim{\rm GPD}_{\sigma_{u},\xi};$
(16)

the distribution function is given in (7). The extremal correlation observed
in Figure 5 is accommodated by allowing $\sigma_{u}$ to depend on POH, MESHS,
exposure covariates and unobserved spatial discrepancies and time-related
autocorrelation, respectively modeled with a Gaussian process and an auto-
regressive process, leading to

$\log\sigma_{u,t}(s)=\sigma_{0}+\sigma_{1}M_{s,t}^{\text{MESHS}}+\sigma_{2}M_{s,t}^{\text{MESHS
$\cdot$ POH}}+\sigma_{3}M_{s}^{\text{Exp}}+X^{\sigma}(s)$ (17)

where $X^{\sigma}$ is a spatial Gaussian process with zero mean and a Matérn
covariance matrix (6), in which the Euclidean distance has been replaced by
the chordal distance because our Gaussian process occurs on the surface of a
sphere, whose curvature should be reflected by our model (Jeong, Jun and
Genton, 2017). The chordal distance is the length of a line passing through
the three-dimensional Earth to connect two points on its surface. For two
locations $s_{1}$ and $s_{2}$ with respective geographical coordinates
$(x_{1},y_{1})$, $(x_{2},y_{2})$, the chordal distance between $s_{1}$ and
$s_{2}$ is defined by

${\rm
C}(s_{1},s_{2})=2r\arcsin\left[\frac{1}{2}\left\\{1-\cos_{p}(y_{2}-y_{1})+\cos_{p}y_{1}\cos_{p}y_{2}(1-\cos_{p}(x_{2}-x_{1})\right\\}\right]^{1/2},$

where $r=6371$km is the Earth’s radius and $\cos_{p}(z)=\cos(\pi z/180)$. In a
small area such as the canton, using the chordal distance instead of the
Euclidean distance might not make a huge difference, but it would matter if
the model was used for larger regions.

In view of Figure 5 we allowed the shape parameter to vary as

$\xi(t)=\begin{cases}\xi_{1},&t\in\\{\text{May, June, July, August}\\},\\\
\xi_{2},&\text{ otherwise.}\end{cases}$ (18)

##  5 Model fitting and validation

###  5.1 Technical challenges

Fitting the Bayesian hierarchical model described in Section 4 is challenging
due to its complexity, the size of the parameter space and the large number of
data points. Recent advances in spatial statistics allow better computational
efficiency for Bayesian models with latent variables (Rue et al., 2017). In
our case, it would be desirable to use R-INLA, which has been widely and
successfully used for environmental data (e.g., Castro-Camilo, Huser and Rue,
2019; Koh et al., 2023). As our work is part of a collaboration involving
several subprojects mostly written in the programming language Python, the
model described in Sections 3 and 4 was also coded in Python so that our
collaborators would find it accessible. There is no equivalent of R-INLA in
Python, and reproducing it for Python users would have taken far too long, so
despite the resulting drop in computational efficiency we resorted to Markov
chain Monte Carlo (MCMC) methods.

In contrast to Metropolis–Hastings steps, which make trajectory proposals
within a possibly skewed ball (Hastings, 1970; Metropolis et al., 2004), or to
Gibbs sampling, which generally only moves in a few dimensions at a time
(Gelfand, 2000), Hamiltonian Monte Carlo (HMC) generates proposals based on
the shape of the posterior by using its gradient (Betancourt, 2017). In MCMC
algorithms, the termination criterion identifies when a trajectory is long
enough for adequate exploration of the neighborhood around the current state,
but in HMC, this criterion should be chosen to compromise between taking full
advantage of the Hamiltonian trajectories and wise use of computational
resources (Betancourt, 2017). The No-U-Turn Sampler (NUTS) is a HMC algorithm
that proves particularly efficient in converging for high-dimensional
posterior distributions (Homan and Gelman, 2014). Indeed, NUTS uses a dynamic
termination criterion that considers only the position and momentum of a
trajectory’s boundaries: when it is met, further sampling typically leads to
neighborhoods that have already been explored. In addition to this specific
termination criterion, NUTS implements a multiplicative expansion of the
trajectory that allows fast exploration of the parameter space within limited
computer memory (Betancourt, 2017). We used NUTS for the count model described
in Section 4.1 and for the extremal model described in Section 4.2.2.

For the non-extreme claims model detailed in Section 4.2.1, a differential
evolution Metropolis (DE-MC) sampling step with a snooker updater was used, as
it is more efficient and faster than the classical random walk Metropolis
step. DE-MC combines a differential evolution genetic algorithm and MCMC
simulation (Ter Braak, 2006). The snooker updater makes it less
computationally expensive than classical DE-MC, as it updates different chains
in parallel with information from past states (Ter Braak and Vrugt, 2008), and
is faster than NUTS for this model, with no significant impact on the results.

We systematically exclude the initial 500 samples drawn, which are reserved
for a tuning phase during which the sampler dynamically adjusts the step sizes
and scalings to optimize its subsequent performance. We monitor the
convergence of the model parameters using informal diagnostic plots; see the
Supplementary Material. We check that autocorrelation has decreased to
approximately zero during the sampling, and examine trace plots of the sampled
parameters for the absence of patterns. Running the claim counts model took
about five hours for about a thousand parameters. For the claim values,
fitting the model took two hours for the GPD model (38 parameters) and less
than an hour for the Beta model (36 parameters). Non-informative priors were
found to perform significantly better than weakly informative priors in our
case and thus were attributed to all of the model parameters. To make sure the
posterior is proper, we check that its distribution percentiles and mean are
finite and reasonable.

###  5.2 Metrics

To assess the model’s performance in improving spatial patterns we use
diagnostic quantities that include the following two specific metrics.

The spatially convolved Kolmogorov–Smirnov statistic (Miralles et al., 2022)
represents the disagreement between the spatial distributions of the generated
and observed images and is computed as the maximum absolute difference of
empirical cumulative distribution functions for the generated and true damage,
summed over $10\times 10$ patches of the image of interest. The aim is to
obtain a metric with properties close to those of the Fréchet inception
distance (Heusel et al., 2017) for images by assessing the match between
predictions and targets, as a human eye would. After extracting $M$ spatial
patches of constant size from the target and predicted images, we set

$\mathrm{SKSS}=\sum_{t\leq N_{T},j\leq
M}\max_{x\in\mathbb{R}}|F_{jt}(x)-\hat{F}_{jt}(x)|,$

where $F_{jt}$ represents the empirical cumulative distribution function of
the hail damage for a single spatial patch $j$ and time $t$ and $\hat{F}_{jt}$
is its analog for predicted damage. This metric evaluates the local agreement
between two distributions rather than focusing on individual pixels.

The log-spectral distance (Rabiner and Juang, 1993) is expressed as the log-
difference of power spectra between the generated and realized samples,

$\mathrm{LSD}=\left\\{\dfrac{1}{2N_{T}\times P}\sum_{t\leq N_{T},i\leq
P}\left[10\log_{10}\left(\dfrac{|g(c_{it})|^{2}}{|g(\hat{c}_{it})|^{2}}\right)\right]^{2}\right\\}^{1/2},$

where $g$ is the Fourier transform, $|g(\cdot)|^{2}$ the power spectrum, $c$
is the target map of damage and $\hat{c}$ its estimate. This evaluates whether
the generated images reproduce the spatial structures noticeable in the target
images.

##  6 Results

There is no visible long-term trend in either the claim count or value in the
insurance data; we can thus split data into sets of consecutive years. The
training set is built from years up to 2015, the validation data comprises the
years 2016–2017, and later years are used as the test set. In the following
analysis of the results, unless specified otherwise, the average prediction
for the test set over 1000 different sets of parameters sampled from the
posterior distribution is used to construct graphs and maps. We recall that
our objectives are to accurately capture spatial patterns for hail damage, to
be able to predict localized and extreme damage and to match the distribution
of the target data provided by GVZ. We shall see that the fitted random line
process with extreme marks achieves this. We start by evaluating the
performance of the Gaussian line process, then explain the procedure for
combining counts and claim values, and finally discuss predicted claims.

###  6.1 Claim counts

Our benchmark for evaluating the performance of the random line model
presented in Section 3.1 is the percentage of affected assets (PAA)-based
gridded claim count predicted with CLIMADA (Schmid et al., 2023+). The PAA,
defined as the per-cell proportion of damaged buildings, is expected to
increase with the value of MESHS, since hailstorms with larger hailstones
should cause more damage. Figure 6(a) shows that the observed PAA does not
increase linearly and is very variable, but that the predicted and observed
values are quite similar, while the 95% prediction range captures the observed
variation well. The impact function computed through CLIMADA tends to over-
predict the percentage of affected assets for any observed MESHS value.

Figure 6(b) suggests that small claim counts are over-predicted by our model.
For days with more than 1000 recorded claims, the distribution of predicted
counts is very close to the observed claim counts; the line process captures
days with very many claims particularly well.

Table 1 assesses how much our model improves on CLIMADA in terms of predicting
the daily claim count. The random line model reduces the false alarm rate by
about 40% and increases the positive predictive value by 15% and the
specificity by 42%, so it makes fewer mistakes on average in predicting both
positive and zero counts. Compared to CLIMADA total predicted counts per day,
Table 1 shows that the sensitivity has dropped by 12%, i.e., our model might
miss days with a positive claim count, but inspection of the data reveals that
it only misses days with fewer than ten claims and less than CHF 10K overall
damage.

Examples of predicted counts plotted in Figure 7 show that the line model
helps to concentrate the predicted damage on straight lines, giving results
that resemble observed claim counts which are usually concentrated in hail
streaks of width just a few km in the Alpine region (Nisi et al., 2018). In
contrast, the predictions from CLIMADA are broadly distributed according to
the MESHS footprint, which typically covers a whole storm cell core (Nisi et
al., 2016). The average predicted count over the canton is also closer to the
realized value using the line model than with CLIMADA.

Figure 6: Comparison of CLIMADA and model: (a) predicted and observed percentages of affected assets; (b) QQ-plots of realized versus predicted quantiles for the number of claims per grid cell per day. Table 1: Comparison of the false alarm rate, sensitivity, specificity and positive predictive value (%) averaged over the time dimension for CLIMADA’s hail damage model and the random line model. If $a$, $b$, $c$ and $d$ denote the true positive, false positive, false negative and true negative numbers, the false alarm rate is computed as $b/(b+d)$, the sensitivity as $a/(a+c)$, the specificity as $d/(b+d)$ and the positive predictive value as $a/(a+b)$. | False Alarm | Sensitivity | Specificity | Positive Predictive Value
---|---|---|---|---
CLIMADA | $72.1$ | $64.8$ | $27.9$ | $62.9$
Model | $29.7$ | $52.1$ | $70.3$ | $77.0$

Figure 7: Comparison of locations of observed claim counts (left), CLIMADA
predicted counts (center), and our predicted counts (right) for three dates
selected over all dates with more than 10 observed claims over the canton on
the 2000–2021 period. The titles give the observed number of claims, CLIMADA
predicted count, and the average count and its 95% predicted range from the
random line model (right).

###  6.2 Combination of claim counts and values

One way to combine counts and values would be to compute a per-cell impact
coefficient corresponding to the predicted number of buildings impacted by a
hailstorm divided by the total number of buildings in the grid cell. This
would then be multiplied by the total possible hail damage in the cell, i.e.,
the sum of predicted values of claims for all buildings, to obtain the per-
cell effective damage. This possibility was considered but not pursued, since
our aim is to predict claim values for single buildings and not the aggregated
damage in a cell. Our combination of counts and values thus involves choosing
which buildings are impacted by hail given the predicted count in a cell. In
each cell, buildings are first sorted by their exposure (i.e., insured value),
and the first $N$ are selected to compute the total damage. To predict damage
and its confidence range, we sample the counts for the cell $n$ times and the
claim value for every building $m$ times, apply the procedure described in the
previous sentence to the $mn$ samples, and finally compute the average total
damage and its 95$\%$ prediction range.

###  6.3 Claim values

Evaluation of the full hail damage prediction involves the combination of
counts and claim values, as described in Section 6.2. We use the two metrics
described in Section 5.2 to compare the spatial patterns of hail damage
predicted with CLIMADA and our model. Figure 8(a) shows that both take higher
values for CLIMADA than for our model, so the latter better captures the
observed hail damage footprints. Figure 8 also shows that the distribution of
hail damage predicted with our model is close to the observed distribution
both on a single building scale (b) and on a 2km square grid (c), though non-
extreme damage is slightly over-predicted, which suggests further research on
the distribution of non-extreme claim values might be needed. CLIMADA
systematically under-predicts the aggregated damage per grid cell, and there
is a clear improvement using the random line process. Figure 8(b) compares the
model’s input, downscaled CLIMADA, to the prediction. As expected, downscaled
CLIMADA under-predicts damage per building (see Section 2.2.2), while the
random line model predicts realistic values, particularly so for claims above
CHF 5000.

Figure 9 shows some daily hail damage maps for days on which there was over
one million Swiss francs of realized hail impacts in the canton (some of these
claim dates belong to the train or validation set). The predicted claim values
appear to be locally large, matching the spatial pattern of realized damage,
whereas CLIMADA damage is more dispersed. The average total predicted damage
over the canton for days with extreme realized hail damage is close to the
observed value, which the confidence interval usually captures. The line model
is thus able to predict well-located extremes while providing reliable
estimates of the total damage. The lowest panel in Figure 9 shows the most
extreme hail event in the two decades of our data, on 28 June 2021, which
involved roof-penetrating hail damage with 177 claims above CHF 100,000. Our
model manages to capture extreme damage on this day, with average predicted
values up to CHF 180,000. The spatial pattern using both CLIMADA and, to a
lesser extent, our model, is wider than the observed data, which might be
related to overestimation of the MESHS intensity that day in the northern half
of the canton (Figure 1).

Figure 8: Comparison of metrics of predicted damage. (a) values of scaled LSD
and SKSS for our model, with those from CLIMADA subtracted. (b) QQ-plots of
realized versus predicted quantiles for the damage per building and (c) per
2km grid cell, with dashed blue lines showing the 95% prediction range.

Figure 9: Example daily hail impact maps. The columns represent the observed
claims (left), CLIMADA-predicted claims (center), and the hail damage
prediction using the random line model (right). The left color bar relates to
the prediction at the scale of the cell (i.e. relevant for CLIMADA predicted
damage), while the right color bar displays a log scale for the per-building
claim values (i.e. relevant for realized and predicted damage). The observed
monetary cantonal damage (left), CLIMADA predicted total (middle), and average
damage from the random line model with its 95% predicted range (right) are
displayed in the titles for each selected date.

## Conclusion

The model developed in this paper seems to be the first to combine a Gaussian
random line process with extreme-value modeling in order to predict the
spatial footprint of hail damage. It improves on the use of a benchmark
deterministic hail damage function: in particular, it captures extreme damage
values for individual buildings well, reproduces the spatial pattern of hail
in the insurance data, and its stochasticity enables uncertainty
quantification. With appropriate changes, such as for instance the possibility
of modeling multiple random lines at the same time, our approach could be
generalized to larger areas and would be useful in studying the insurance
impacts of climate change. It would be interesting to use thunderstorm cell
direction instead of large-scale wind direction as the covariate for the slope
of the random line process.

## Authorship contribution statement

##### Ophélia Miralles

Methodology, Software, Validation, Formal Analysis, Investigation, Data
Curation, Writing (original draft, review, editing), Visualization.

##### Anthony Davison

Methodology, Supervision, Editing, Funding acquisition.

##### Timo Schmid

Data Curation, Editing.

[Acknowledgments] We thank the Federal Office of Meteorology and Climatology
(MeteoSwiss) for providing the POH and MESHS data, the insurance company GVZ
for providing the building damage and exposure data, and Daniel Steinfeld of
GVZ for valuable insights.

Ophélia Miralles’ work was partially funded by the Swiss National Science
Foundation (project 200021_178824). Timo Schmid acknowledges funding from the
Swiss National Science Foundation (SNSF) Sinergia grant (CRSII5_201792).

## Data statement

The exploratory analysis, model, and diagnostics are implemented in Python and
the code is freely available on GitHub (github.com/OpheliaMiralles/hail-
damage-modeling). MESHS and POH data are available from MeteoSwiss on demand.
The GVZ insurance data are private, and as such are available for use only
within the scClim project for research purposes. The CLIMADA impact function
is open-source and can be run using the GitHub repository github.com/CLIMADA-
project/climada_papers. ERA5 reanalysis data can be downloaded freely from the
Copernicus Climate Data Store (climate.copernicus.eu/climate-reanalysis).

## References

  * Aznar-Siguan and Bresch (2019) [author] Aznar-Siguan, G.G. and Bresch, D. N.D. N. (2019). CLIMADA v1: a global weather and climate risk assessment platform. Geoscientific Model Development 12 3085–3097. 10.5194/gmd-12-3085-2019
  * Betancourt (2017) [author] Betancourt, MichaelM. (2017). A conceptual introduction to Hamiltonian Monte Carlo. arXiv e-prints arXiv:1701.02434. 10.48550/arXiv.1701.02434
  * Botzen, Bouwer and van den Bergh (2010) [author] Botzen, W. J. W.W. J. W., Bouwer, L. M.L. M. and van den Bergh, J. C. J. M.J. C. J. M. (2010). Climate change and hailstorm damage: Empirical evidence and implications for agriculture and insurance. Resource and Energy Economics 32 341–362. https://doi.org/10.1016/j.reseneeco.2009.10.004
  * Brown, Pogorzelski and Giammanco (2015) [author] Brown, Tanya MT. M., Pogorzelski, William HW. H. and Giammanco, Ian MI. M. (2015). Evaluating hail damage using property insurance claims data. Weather, Climate, and Society 7 197–210.
  * Castro-Camilo, Huser and Rue (2019) [author] Castro-Camilo, DanielaD., Huser, RaphaëlR. and Rue, HåvardH. (2019). A spliced gamma-generalized Pareto model for short-term extreme wind speed probabilistic forecasting. Journal of Agricultural, Biological and Environmental Statistics 24 517–534.
  * Changnon (2008) [author] Changnon, Stanley A.S. A. (2008). Temporal and spatial distributions of damaging hail in the continental United States. Physical Geography 29 341–350. 10.2747/0272-3646.29.4.341
  * Deepen (2006) [author] Deepen, JanJ. (2006). Schadenmodellierung extremer Hagelereignisse in Deutschland, Master’s thesis, Westfälischen Wilhelms-Universität Münster.
  * Gelfand (2000) [author] Gelfand, Alan E.A. E. (2000). Gibbs sampling. Journal of the American Statistical Association 95 1300–1304. 10.1080/01621459.2000.10474335
  * Hastings (1970) [author] Hastings, W. K.W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57 97–109. 10.1093/biomet/57.1.97
  * Hersbach et al. (2020) [author] Hersbach, HansH., Bell, BillB., Berrisford, PaulP., Hirahara, ShojiS., Horányi, AndrásA., Muñoz-Sabater, JoaquínJ., Nicolas, JulienJ., Peubey, CaroleC., Radu, RalucaR., Schepers, DinandD. et al. (2020). The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society 146 1999–2049.
  * Heusel et al. (2017) [author] Heusel, MartinM., Ramsauer, HubertH., Unterthiner, ThomasT., Nessler, BernhardB. and Hochreiter, SeppS. (2017). GANs trained by a two time-Scale update rule converge to a local Nash equilibrium. In Advances in Neural Information Processing Systems (I.I. Guyon, U. VonU. V. Luxburg, S.S. Bengio, H.H. Wallach, R.R. Fergus, S.S. Vishwanathan and R.R. Garnett, eds.) 30 6626–6637. Curran Associates, Inc., New Orleans, Louisiana, United States.
  * Hohl, Schiesser and Aller (2002) [author] Hohl, RomanR., Schiesser, Hans-HeinrichH.-H. and Aller, DörteD. (2002). Hailfall: the relationship between radar-derived hail kinetic energy and hail damage to buildings. Atmospheric Research 63 177–207. https://doi.org/10.1016/S0169-8095(02)00059-5
  * Homan and Gelman (2014) [author] Homan, Matthew D.M. D. and Gelman, AndrewA. (2014). The No-U-Turn Sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research 15 1593–1623. 3214779
  * Jeong, Jun and Genton (2017) [author] Jeong, JaehongJ., Jun, MikyoungM. and Genton, Marc GM. G. (2017). Spherical process models for global spatial statistics. Statistical Science 32 501–513.
  * Koh et al. (2023) [author] Koh, JonathanJ., Pimont, FrançoisF., Dupuy, Jean-LucJ.-L. and Opitz, ThomasT. (2023). Spatiotemporal wildfire modeling through point processes with moderate and extreme marks. Annals of Applied Statistics 17 560–582. 10.1214/22-AOAS1642
  * Kopp et al. (2023) [author] Kopp, JérômeJ., Schröer, KatharinaK., Schwierz, CorneliaC., Hering, AlessandroA., Germann, UrsU. and Martius, OliviaO. (2023). The summer 2021 Switzerland hailstorms: weather situation, major impacts and unique observational data. Weather 78 184–191. https://doi.org/10.1002/wea.4306
  * Liu, Li and Wang (2021) [author] Liu, KaiK., Li, PeiP. and Wang, ZhonggangZ. (2021). Statistical modeling of random hail impact. Extreme Mechanics Letters 48 101374\. https://doi.org/10.1016/j.eml.2021.101374
  * Metropolis et al. (2004) [author] Metropolis, NicholasN., Rosenbluth, Arianna W.A. W., Rosenbluth, Marshall N.M. N., Teller, Augusta H.A. H. and Teller, EdwardE. (2004). Equation of state calculations by fast computing machines. The Journal of Chemical Physics 21 1087–1092. 10.1063/1.1699114
  * Miralles et al. (2022) [author] Miralles, OphéliaO., Steinfeld, DanielD., Martius, OliviaO. and Davison, Anthony CA. C. (2022). Downscaling of historical wind fields over Switzerland using generative adversarial networks. Artificial Intelligence for the Earth Systems 1 e220018.
  * Müller (2021) [author] Müller, AndréA. (2021). 2 Milliarden Franken wegen Sturm und Hagel: 2021 wird für die Schweiz eines der teuersten Schadenjahre aller Zeiten. https://www.nzz.ch/wirtschaft/unwetter-schweiz-2021-hagel-und-sturm-kosteten-2-mrd-franken-ld.1652483?reduced=true. Accessed:2023-06-25.
  * Nisi et al. (2016) [author] Nisi, LucaL., Martius, OliviaO., Hering, AlessandroA., Kunz, MM. and Germann, UU. (2016). Spatial and temporal distribution of hailstorms in the Alpine region: a long-term, high resolution, radar-based analysis. Quarterly Journal of the Royal Meteorological Society 142 1590–1604.
  * Nisi et al. (2018) [author] Nisi, LucaL., Hering, AlessandroA., Germann, UrsU. and Martius, OliviaO. (2018). A 15-year hail streak climatology for the Alpine region. Quarterly Journal of the Royal Meteorological Society 144 1429-1449. https://doi.org/10.1002/qj.3286
  * Otto (2009) [author] Otto, MM. (2009). Modellierung von Hagelschäden in der Pkw-Kaskoversicherung in Deutschland, Master’s thesis, Technische Universität Dresden.
  * Perera et al. (2018) [author] Perera, ShiharaS., Lam, NelsonN., Pathirana, MahilM., Zhang, LihaiL., Ruan, DongD. and Gad, EmadE. (2018). Probabilistic modelling of forces of hail. Natural Hazards 91 133–153. 10.1007/s11069-017-3117-7
  * Punge et al. (2014) [author] Punge, HJH., Bedka, KMK., Kunz, MM. and Werner, AA. (2014). A new physically based stochastic event catalog for hail in Europe. Natural Hazards 73 1625–1645.
  * Púčik et al. (2017) [author] Púčik, TomášT., Groenemeijer, PieterP., Rädler, Anja T.A. T., Tijssen, LarsL., Nikulin, GrigoryG., Prein, Andreas F.A. F., van Meijgaard, ErikE., Fealy, RowanR., Jacob, DanielaD. and Teichmann, ClaasC. (2017). Future changes in European severe convection environments in a regional climate model ensemble. Journal of Climate 30 6771–6794. https://doi.org/10.1175/JCLI-D-16-0777.1
  * Rabiner and Juang (1993) [author] Rabiner, L. R.L. R. and Juang, B. H.B. H. (1993). Fundamentals of Speech Recognition. Prentice Hall PTR, Englewood Cliffs, New Jersey.
  * Rädler et al. (2019) [author] Rädler, Anja T.A. T., Groenemeijer, Pieter H.P. H., Faust, EberhardE., Sausen, RobertR. and Púčik, TomášT. (2019). Frequency of severe thunderstorms across Europe expected to increase in the 21st century due to rising instability. npj Climate and Atmospheric Science 2 30\. 10.1038/s41612-019-0083-7
  * Rue et al. (2017) [author] Rue, HåvardH., Riebler, AndreaA., Sørbye, Sigrunn HS. H., Illian, Janine BJ. B., Simpson, Daniel PD. P. and Lindgren, Finn KF. K. (2017). Bayesian computing with INLA: a review. Annual Review of Statistics and Its Application 4 395–421.
  * Schmid et al. (2023+) [author] Schmid, TimoT., Villiger, LeonieL., Portmann, RaphaelR. and Bresch, David N.D. N. (2023+). Open-source hail damage model for buildings and cars. Natural Hazards and Earth System Sciences.
  * Schmidberger, Daniell and Kunz (2018) [author] Schmidberger, MM., Daniell, JEJ. and Kunz, MM. (2018). Hail hazard and hail risk modeling for Germany based on a combination of radar and insurance data. In 1st North American Workshop on Hail & Hailstorms, 14–16 August 2018 (Andrew J.A. J. Heymsfield and Ian M.I. M. Giammanco, eds.). American Meteorological Society, Boulder, Colorado.
  * Schuster, Blong and McAneney (2006) [author] Schuster, Sandra S.S. S., Blong, Russell J.R. J. and McAneney, K. JohnK. J. (2006). Relationship between radar-derived hail kinetic energy and damage to insured buildings for severe hailstorms in Eastern Australia. Atmospheric Research 81 215–235. https://doi.org/10.1016/j.atmosres.2005.12.003
  * Ter Braak (2006) [author] Ter Braak, Cajo JFC. J. (2006). A Markov chain Monte Carlo version of the genetic algorithm differential evolution: easy Bayesian computing for real parameter spaces. Statistics and Computing 16 239–249.
  * Ter Braak and Vrugt (2008) [author] Ter Braak, Cajo JFC. J. and Vrugt, Jasper AJ. A. (2008). Differential evolution Markov chain with snooker updater and fewer chains. Statistics and Computing 18 435–446.
  * Varty et al. (2021) [author] Varty, ZakZ., Tawn, Jonathan AJ. A., Atkinson, Peter MP. M. and Bierman, StijnS. (2021). Inference for extreme earthquake magnitudes accounting for a time-varying measurement process. arXiv preprint arXiv:2102.00884.
  * Warren et al. (2020) [author] Warren, Robert A.R. A., Ramsay, Hamish A.H. A., Siems, Steven T.S. T., Manton, Michael J.M. J., Peter, Justin R.J. R., Protat, AlainA. and Pillalamarri, AnuA. (2020). Radar-based climatology of damaging hailstorms in Brisbane and Sydney, Australia. Quarterly Journal of the Royal Meteorological Society 146 505–530. https://doi.org/10.1002/qj.3693

Additional details Here we provide additional information and figures about
model selection, including material about the choice of fixed hyperparameters
such as the threshold for the GPD model, and also MC diagnostics related to
the validation of the Bayesian model parameters.

Figure 10: Threshold selection method described in Varty et al. (2021) applied
to the log total damage per cell. The left panel shows a minimum
qq-$\ell_{1}$-distance for a threshold of 8.06. The QQ-plot for the GPD fit of
the log total sum of damage above this threshold is displayed in the right
panel, in which the profile likelihood-based 95% confidence interval is shown
by the blue shaded area. Figure 11: Evolution of the autocorrelation through
sampling for parameters of the Negative Binomial model presented in Section
4.1. The grey area designates the acceptable range for the autocorrelation at
the end of sampling to assume convergence of the model. The bounds of the
confidence range are computed from the central limit theorem. Figure 12:
Autocorrelation plot (left) and trace plot (right) for the posterior
distribution of the shape parameter $\alpha$ in the Negative Binomial model
presented in Section 4.1. A close-to-zero autocorrelation through sampling and
no specific trend or pattern in the trace plot is usually a good sign of model
convergence. Figure 13: Kernel density estimate plot: (left) and trace plot:
(right) for parameters of the Beta model presented in Section 4.2.1. No
specific trend or pattern in the trace plot is usually a good sign of model
convergence. Figure 14: Evolution of the autocorrelation through sampling for
parameters of the GPD model presented in Section 4.2.2. The grey area
designates the acceptable range for the autocorrelation at the end of sampling
to assume convergence of the model. The bounds of the confidence range are
computed from the central limit theorem. Figure 15: Autocorrelation plot
(left) and trace plot (right) for the posterior distribution of the shape
parameter $\xi$ in the GPD model presented in Section 4.2.2. A close-to-zero
autocorrelation through sampling and no specific trend or pattern in the trace
plot is usually a good sign of model convergence.


* (50) T. Modak, J. C. Romão, S. Sadhukhan, J. a. P. Silva, and R. Srivastava, Constraining wrong-sign $hbb$ couplings with $h\rightarrow\Upsilon\gamma$, Phys. Rev. D 94 (2016), no. 7 075017, [arXiv:1607.07876].
  * (51) M. Raju, J. P. Saha, D. Das, and A. Kundu, Double Higgs boson production as an exclusive probe for a sequential fourth generation with wrong-sign Yukawa couplings, Phys. Rev. D 101 (2020), no. 5 055036, [arXiv:2001.05280].
  * (52) F. U. Bernlochner, Z. Ligeti, M. Papucci, and D. J. Robinson, Combined analysis of semileptonic $B$ decays to $D$ and $D^{*}$: $R(D^{(*)})$, $|V_{cb}|$, and new physics, Phys. Rev. D 95 (2017), no. 11 115008, [arXiv:1703.05330]. [Erratum: Phys.Rev.D 97, 059902 (2018)].
  * (53) I. Caprini, L. Lellouch, and M. Neubert, Dispersive bounds on the shape of anti-B —$>$ D(*) lepton anti-neutrino form-factors, Nucl. Phys. B 530 (1998) 153–181, [hep-ph/9712417].
  * (54) Y. Sakaki, M. Tanaka, A. Tayduganov, and R. Watanabe, Testing leptoquark models in $\bar{B}\to D^{(*)}\tau\bar{\nu}$, Phys. Rev. D 88 (2013), no. 9 094012, [arXiv:1309.0301].
  * (55) A. Bharucha, D. M. Straub, and R. Zwicky, $B\to V\ell^{+}\ell^{-}$ in the Standard Model from light-cone sum rules, JHEP 08 (2016) 098, [arXiv:1503.05534].
  * (56) N. Gubernari, A. Kokulu, and D. van Dyk, $B\to P$ and $B\to V$ Form Factors from $B$-Meson Light-Cone Sum Rules beyond Leading Twist, JHEP 01 (2019) 150, [arXiv:1811.00983].
  * (57) W. Detmold and S. Meinel, $\Lambda_{b}\to\Lambda\ell^{+}\ell^{-}$ form factors, differential branching fraction, and angular observables from lattice QCD with relativistic $b$ quarks, Phys. Rev. D 93 (2016), no. 7 074501, [arXiv:1602.01399].
  * (58) V. Bernard, M. Oertel, E. Passemar, and J. Stern, K(mu3)**L decay: A Stringent test of right-handed quark currents, Phys. Lett. B 638 (2006) 480–486, [hep-ph/0603202].
  * (59) V. Bernard, M. Oertel, E. Passemar, and J. Stern, Dispersive representation and shape of the K(l3) form factors: Robustness, Phys. Rev. D 80 (2009) 034034, [arXiv:0903.1654].
  * (60) FlaviaNet Working Group on Kaon Decays Collaboration, M. Antonelli et al., An Evaluation of $|V_{us}|$ and precise tests of the Standard Model from world data on leptonic and semileptonic kaon decays, Eur. Phys. J. C 69 (2010) 399–424, [arXiv:1005.2323].
  * (61) LHCb Collaboration, R. Aaij et al., Measurement of $|V_{cb}|$ with $B_{s}^{0}\to D_{s}^{(*)-}\mu^{+}\nu_{\mu}$ decays, Phys. Rev. D 101 (2020), no. 7 072004, [arXiv:2001.03225].
  * (62) E. McLean, C. Davies, J. Koponen, and A. Lytle, $B_{s}\to D_{s}\ell\nu$ Form Factors for the full $q^{2}$ range from Lattice QCD with non-perturbatively normalized currents, Phys. Rev. D 101 (2020), no. 7 074513, [arXiv:1906.00701].
  * (63) HFLAV Collaboration, Y. S. Amhis et al., Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of 2018, arXiv:1909.12524.
  * (64) L. Di Luzio, M. Kirk, A. Lenz, and T. Rauh, $\Delta M_{s}$ theory precision confronts flavour anomalies, JHEP 12 (2019) 009, [arXiv:1909.11087].
  * (65) A. Crivellin, D. Müller, and C. Wiegand, $b\to s\ell^{+}\ell^{-}$ transitions in two-Higgs-doublet models, JHEP 06 (2019) 119, [arXiv:1903.10440].
  * (66) C. Geng and J. N. Ng, Charged Higgs Effect in $B$)$d$)0 - $\bar{B}(d$)0 Mixing, $K\to\pi$ Neutrino Anti-neutrino Decay and Rare Decays of $B$ Mesons, Phys. Rev. D 38 (1988) 2857. [Erratum: Phys.Rev.D 41, 1715 (1990)].
  * (67) J. Urban, F. Krauss, U. Jentschura, and G. Soff, Next-to-leading order QCD corrections for the B0 anti-B0 mixing with an extended Higgs sector, Nucl. Phys. B 523 (1998) 40–58, [hep-ph/9710245].
  * (68) D. King, A. Lenz, and T. Rauh, Bs mixing observables and |Vtd/Vts| from sum rules, JHEP 05 (2019) 034, [arXiv:1904.00940].
  * (69) M. Kirk, A. Lenz, and T. Rauh, Dimension-six matrix elements for meson mixing and lifetimes from sum rules, JHEP 12 (2017) 068, [arXiv:1711.02100]. [Erratum: JHEP 06, 162 (2020)].
  * (70) A. G. Grozin, R. Klein, T. Mannel, and A. A. Pivovarov, $B^{0}-\bar{B}^{0}$ mixing at next-to-leading order, Phys. Rev. D 94 (2016), no. 3 034024, [arXiv:1606.06054].
  * (71) R. J. Dowdall, C. T. H. Davies, R. R. Horgan, G. P. Lepage, C. J. Monahan, J. Shigemitsu, and M. Wingate, Neutral B-meson mixing from full lattice QCD at the physical point, Phys. Rev. D 100 (2019), no. 9 094508, [arXiv:1907.01025].
  * (72) RBC/UKQCD Collaboration, P. A. Boyle, L. Del Debbio, N. Garron, A. Juttner, A. Soni, J. T. Tsang, and O. Witzel, SU(3)-breaking ratios for $D_{(s)}$ and $B_{(s)}$ mesons, arXiv:1812.08791.
  * (73) Fermilab Lattice, MILC Collaboration, A. Bazavov et al., $B^{0}_{(s)}$-mixing matrix elements from lattice QCD for the Standard Model and beyond, Phys. Rev. D 93 (2016), no. 11 113016, [arXiv:1602.03560].
  * (74) A. J. Buras, M. Jamin, and P. H. Weisz, Leading and Next-to-leading QCD Corrections to $\epsilon$ Parameter and $B^{0}-\bar{B}^{0}$ Mixing in the Presence of a Heavy Top Quark, Nucl. Phys. B 347 (1990) 491–536.
  * (75) LHCb Collaboration, R. Aaij et al., Precise determination of the $B^{0}_{s}$-$\overline{B}^{0}_{s}$ oscillation frequency, arXiv:2104.04421.
  * (76) M. Misiak, A. Rehman, and M. Steinhauser, Towards $\overline{B}\to{X}_{s}\gamma$ at the NNLO in QCD without interpolation in mc, JHEP 06 (2020) 175, [arXiv:2002.01548].
  * (77) M. Misiak and M. Steinhauser, NNLO QCD corrections to the anti-B —$>$ X(s) gamma matrix elements using interpolation in m(c), Nucl. Phys. B 764 (2007) 62–82, [hep-ph/0609241].
  * (78) M. Misiak et al., Updated NNLO QCD predictions for the weak radiative B-meson decays, Phys. Rev. Lett. 114 (2015), no. 22 221801, [arXiv:1503.01789].
  * (79) CLEO Collaboration, S. Chen et al., Branching fraction and photon energy spectrum for $b\to s\gamma$, Phys. Rev. Lett. 87 (2001) 251807, [hep-ex/0108032].
  * (80) BaBar Collaboration, J. P. Lees et al., Precision Measurement of the $B\to X_{s}\gamma$ Photon Energy Spectrum, Branching Fraction, and Direct CP Asymmetry $A_{CP}(B\to X_{s+d}\gamma)$, Phys. Rev. Lett. 109 (2012) 191801, [arXiv:1207.2690].
  * (81) Belle Collaboration, A. Abdesselam et al., Measurement of the inclusive $B\to X_{s+d}\gamma$ branching fraction, photon energy spectrum and HQE parameters, in 38th International Conference on High Energy Physics, 8, 2016.  arXiv:1608.02344.
  * (82) F. Borzumati and C. Greub, Two Higgs doublet model predictions for anti-B —$>$ X(s) gamma in NLO QCD: Addendum, Phys. Rev. D 59 (1999) 057501, [hep-ph/9809438].
  * (83) T. Hermann, M. Misiak, and M. Steinhauser, $\bar{B}\to X_{s}\gamma$ in the Two Higgs Doublet Model up to Next-to-Next-to-Leading Order in QCD, JHEP 11 (2012) 036, [arXiv:1208.2788].
  * (84) W. Altmannshofer and P. Stangl, New Physics in Rare B Decays after Moriond 2021, arXiv:2103.13370.
  * (85) LHCb Collaboration, R. Aaij et al., Measurement of the $B^{0}_{s}\to\mu^{+}\mu^{-}$ branching fraction and effective lifetime and search for $B^{0}\to\mu^{+}\mu^{-}$ decays, Phys. Rev. Lett. 118 (2017), no. 19 191801, [arXiv:1703.05747].
  * (86) ATLAS Collaboration, M. Aaboud et al., Study of the rare decays of $B^{0}_{s}$ and $B^{0}$ mesons into muon pairs using data collected during 2015 and 2016 with the ATLAS detector, JHEP 04 (2019) 098, [arXiv:1812.03017].
  * (87) CMS Collaboration, A. M. Sirunyan et al., Measurement of properties of B${}^{0}_{\mathrm{s}}\to\mu^{+}\mu^{-}$ decays and search for B${}^{0}\to\mu^{+}\mu^{-}$ with the CMS experiment, JHEP 04 (2020) 188, [arXiv:1910.12127].
  * (88) LHCb Collaboration Collaboration, Combination of the ATLAS, CMS and LHCb results on the $B^{0}_{(s)}\to\mu^{+}\mu^{-}$ decays, Tech. Rep. LHCb-CONF-2020-002. CERN-LHCb-CONF-2020-002, CERN, Geneva, Aug, 2020.
  * (89) LHCb Collaboration, R. Aaij et al., Measurement of the $B^{0}_{s}\to\mu^{+}\mu^{-}$ decay properties and search for the $B^{0}\to\mu^{+}\mu^{-}$ and $B^{0}_{s}\to\mu^{+}\mu^{-}\gamma$ decays, arXiv:2108.09283.
  * (90) LHCb Collaboration, R. Aaij et al., Analysis of neutral $B$-meson decays into two muons, arXiv:2108.09284.
  * (91) G. Buchalla and A. J. Buras, QCD corrections to rare K and B decays for arbitrary top quark mass, Nucl. Phys. B 400 (1993) 225–239.
  * (92) C. Bobeth, M. Gorbahn, T. Hermann, M. Misiak, E. Stamou, and M. Steinhauser, $B_{s,d}\to l^{+}l^{-}$ in the Standard Model with Reduced Theoretical Uncertainty, Phys. Rev. Lett. 112 (2014) 101801, [arXiv:1311.0903].
  * (93) M. Beneke, C. Bobeth, and R. Szafron, Power-enhanced leading-logarithmic QED corrections to $B_{q}\to\mu^{+}\mu^{-}$, JHEP 10 (2019) 232, [arXiv:1908.07011].
  * (94) ETM Collaboration, A. Bussone et al., Mass of the b quark and B -meson decay constants from Nf=2+1+1 twisted-mass lattice QCD, Phys. Rev. D 93 (2016), no. 11 114505, [arXiv:1603.04306].
  * (95) A. Bazavov et al., $B$\- and $D$-meson leptonic decay constants from four-flavor lattice QCD, Phys. Rev. D 98 (2018), no. 7 074512, [arXiv:1712.09262].
  * (96) C. Hughes, C. T. H. Davies, and C. J. Monahan, New methods for B meson decay constants and form factors from lattice NRQCD, Phys. Rev. D 97 (2018), no. 5 054509, [arXiv:1711.09981].
  * (97) H. E. Logan and U. Nierste, $B_{s,d}\to\ell^{+}\ell^{-}$ in a two Higgs doublet model, Nucl. Phys. B 586 (2000) 39–55, [hep-ph/0004139].
  * (98) C. Bobeth, G. Hiller, and G. Piranishvili, Angular distributions of $\bar{B}\to\bar{K}\ell^{+}\ell^{-}$ decays, JHEP 12 (2007) 040, [arXiv:0709.4174].
  * (99) M. Bordone, G. Isidori, and A. Pattori, On the Standard Model predictions for $R_{K}$ and $R_{K^{*}}$, Eur. Phys. J. C 76 (2016), no. 8 440, [arXiv:1605.07633].
  * (100) LHCb Collaboration, R. Aaij et al., Test of lepton universality in beauty-quark decays, arXiv:2103.11769.
  * (101) A. Khodjamirian, T. Mannel, A. A. Pivovarov, and Y. M. Wang, Charm-loop effect in $B\to K^{(*)}\ell^{+}\ell^{-}$ and $B\to K^{*}\gamma$, JHEP 09 (2010) 089, [arXiv:1006.4945].
  * (102) A. Khodjamirian, T. Mannel, and Y. M. Wang, $B\to K\ell^{+}\ell^{-}$ decay at large hadronic recoil, JHEP 02 (2013) 010, [arXiv:1211.0234].
  * (103) A. Khodjamirian and A. V. Rusov, $B_{s}\to K\ell\nu_{\ell}$ and $B_{(s)}\to\pi(K)\ell^{+}\ell^{-}$ decays at large recoil and CKM matrix elements, JHEP 08 (2017) 112, [arXiv:1703.04765].
  * (104) M. Algueró, B. Capdevila, A. Crivellin, S. Descotes-Genon, P. Masjuan, J. Matias, M. Novoa Brunet, and J. Virto, Emerging patterns of New Physics with and without Lepton Flavour Universal contributions, Eur. Phys. J. C 79 (2019), no. 8 714, [arXiv:1903.09578]. [Addendum: Eur.Phys.J.C 80, 511 (2020)].
  * (105) T. Hurth, F. Mahmoudi, and S. Neshatpour, Model independent analysis of the angular observables in $B^{0}\to K^{*0}\mu^{+}\mu^{-}$ and $B^{+}\to K^{*+}\mu^{+}\mu^{-}$, Phys. Rev. D 103 (2021) 095020, [arXiv:2012.12207].
  * (106) A. K. Alok, A. Dighe, S. Gangal, and D. Kumar, Continuing search for new physics in $b\to s\mu\mu$ decays: two operators at a time, JHEP 06 (2019) 089, [arXiv:1903.09617].
  * (107) M. Ciuchini, A. M. Coutinho, M. Fedele, E. Franco, A. Paul, L. Silvestrini, and M. Valli, New Physics in $b\to s\ell^{+}\ell^{-}$ confronts new data on Lepton Universality, Eur. Phys. J. C 79 (2019), no. 8 719, [arXiv:1903.09632].
  * (108) T. Hurth, F. Mahmoudi, and S. Neshatpour, Implications of the new LHCb angular analysis of $B\to K^{*}\mu^{+}\mu^{-}$ : Hadronic effects or new physics?, Phys. Rev. D 102 (2020), no. 5 055001, [arXiv:2006.04213].
  * (109) T. Hurth, F. Mahmoudi, D. M. Santos, and S. Neshatpour, More Indications for Lepton Nonuniversality in $b\to s\ell^{+}\ell^{-}$, arXiv:2104.10058.
  * (110) M. Algueró, B. Capdevila, S. Descotes-Genon, J. Matias, and M. Novoa-Brunet, $\boldsymbol{b\to s\ell\ell}$ global fits after Moriond 2021 results, in 55th Rencontres de Moriond on QCD and High Energy Interactions, 4, 2021.  arXiv:2104.08921.
  * (111) C. Cornella, D. A. Faroughy, J. Fuentes-Martín, G. Isidori, and M. Neubert, Reading the footprints of the B-meson flavor anomalies, arXiv:2103.16558.
  * (112) L.-S. Geng, B. Grinstein, S. Jäger, S.-Y. Li, J. Martin Camalich, and R.-X. Shi, Implications of new evidence for lepton-universality violation in $b\to s\ell^{+}\ell^{-}$ decays, arXiv:2103.12738.
  * (113) M. Ciuchini, M. Fedele, E. Franco, A. Paul, L. Silvestrini, and M. Valli, Lessons from the $B^{0,+}\to K^{*0,+}\mu^{+}\mu^{-}$ angular analyses, Phys. Rev. D 103 (2021), no. 1 015030, [arXiv:2011.01212].
  * (114) F. Munir Bhutta, Z.-R. Huang, C.-D. Lü, M. A. Paracha, and W. Wang, New Physics in $b\to s\ell\ell$ anomalies and its implications for the complementary neutral current decays, arXiv:2009.03588.
  * (115) A. Biswas, S. Nandi, I. Ray, and S. K. Patra, New physics in $b\to s\ell\ell$ decays with complex Wilson coefficients, arXiv:2004.14687.
  * (116) P. Athron, C. Balázs, D. H. Jacob, W. Kotlarski, D. Stöckinger, and H. Stöckinger-Kim, New physics explanations of $a_{\mu}$ in light of the FNAL muon $g-2$ measurement, arXiv:2104.03691.
  * (117) F. J. Botella, F. Cornet-Gomez, and M. Nebot, Electron and muon $g-2$ anomalies in general flavour conserving two Higgs doublets models, Phys. Rev. D 102 (2020), no. 3 035023, [arXiv:2006.01934].
  * (118) S. Jana, V. P. K., and S. Saad, Resolving electron and muon $g-2$ within the 2HDM, Phys. Rev. D 101 (2020), no. 11 115037, [arXiv:2003.03386].
  * (119) N. Ghosh and J. Lahiri, Revisiting Generalized Two Higgs Doublet Model in the Light of Muon Anomaly and Lepton Flavor Violating Decays at HL-LHC, arXiv:2010.03590.
  * (120) V. Ilisie, New Barr-Zee contributions to $\mathbf{(g-2)_{\mu}}$ in two-Higgs-doublet models, JHEP 04 (2015) 077, [arXiv:1502.04199].
  * (121) S. M. Barr and A. Zee, Electric Dipole Moment of the Electron and of the Neutron, Phys. Rev. Lett. 65 (1990) 21–24. [Erratum: Phys.Rev.Lett. 65, 2920 (1990)].
  * (122) K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov, Is there a hot electroweak phase transition at m(H) larger or equal to m(W)?, Phys. Rev. Lett. 77 (1996) 2887–2890, [hep-ph/9605288].
  * (123) W. Su, A. G. Williams, and M. Zhang, Strong first order electroweak phase transition in 2HDM confronting future Z & Higgs factories, arXiv:2011.04540.
  * (124) G. C. Dorsch, S. J. Huber, T. Konstandin, and J. M. No, A Second Higgs Doublet in the Early Universe: Baryogenesis and Gravitational Waves, JCAP 05 (2017) 052, [arXiv:1611.05874].
  * (125) L. Wang, Inflation, electroweak phase transition, and Higgs searches at the LHC in the two-Higgs-doublet model, arXiv:2105.02143.
  * (126) P. S. Bhupal Dev and A. Pilaftsis, Maximally Symmetric Two Higgs Doublet Model with Natural Standard Model Alignment, JHEP 12 (2014) 024, [arXiv:1408.3405]. [Erratum: JHEP 11, 147 (2015)].
  * (127) CMS Collaboration, A. M. Sirunyan et al., Measurements of Higgs boson production cross sections and couplings in the diphoton decay channel at $\sqrt{s}=$ 13 TeV, arXiv:2103.06956.
  * (128) ATLAS Collaboration, Measurements of gluon fusion and vector-boson-fusion production of the Higgs boson in $H\rightarrow WW^{*}\rightarrow e\nu\mu\nu$ decays using $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, Tech. Rep. ATLAS-CONF-2021-014, CERN, Geneva, Mar, 2021.
  * (129) ATLAS Collaboration, G. Aad et al., Measurements of Higgs bosons decaying to bottom quarks from vector boson fusion production with the ATLAS experiment at $\sqrt{s}=13\,\text{TeV}$, Eur. Phys. J. C 81 (2021), no. 6 537, [arXiv:2011.08280].
  * (130) ATLAS Collaboration, G. Aad et al., Measurements of $WH$ and $ZH$ production in the $H\rightarrow b\bar{b}$ decay channel in $pp$ collisions at 13 TeV with the ATLAS detector, Eur. Phys. J. C 81 (2021), no. 2 178, [arXiv:2007.02873].
  * (131) LHCb Collaboration, R. Aaij et al., Precise measurement of the $f_{s}/f_{d}$ ratio of fragmentation fractions and of $B^{0}_{s}$ decay branching fractions, arXiv:2103.06810.
  * (132) LHCb Collaboration, R. Aaij et al., Differential branching fractions and isospin asymmetries of $B\to K^{(*)}\mu^{+}\mu^{-}$ decays, JHEP 06 (2014) 133, [arXiv:1403.8044].
  * (133) LHCb Collaboration, R. Aaij et al., Angular analysis of the $B^{+}\rightarrow K^{\ast+}\mu^{+}\mu^{-}$ decay, arXiv:2012.13241.
  * (134) LHCb Collaboration, R. Aaij et al., Measurement of $CP$-Averaged Observables in the $B^{0}\rightarrow K^{*0}\mu^{+}\mu^{-}$ Decay, Phys. Rev. Lett. 125 (2020), no. 1 011802, [arXiv:2003.04831].
  * (135) ATLAS Collaboration, M. Aaboud et al., Angular analysis of $B^{0}_{d}\rightarrow K^{*}\mu^{+}\mu^{-}$ decays in $pp$ collisions at $\sqrt{s}=8$ TeV with the ATLAS detector, JHEP 10 (2018) 047, [arXiv:1805.04000].
  * (136) CMS Collaboration, V. Khachatryan et al., Angular analysis of the decay $B^{0}\to K^{*0}\mu^{+}\mu^{-}$ from pp collisions at $\sqrt{s}=8$ TeV, Phys. Lett. B 753 (2016) 424–448, [arXiv:1507.08126].
  * (137) CMS Collaboration, A. M. Sirunyan et al., Measurement of angular parameters from the decay $\mathrm{B}^{0}\to\mathrm{K}^{*0}\mu^{+}\mu^{-}$ in proton-proton collisions at $\sqrt{s}=$ 8 TeV, Phys. Lett. B 781 (2018) 517–541, [arXiv:1710.02846].
  * (138) Belle Collaboration, A. Abdesselam et al., Test of lepton flavor universality in $B\to K\ell^{+}\ell^{-}$ decays, arXiv:1908.01848.
  * (139) Belle Collaboration, A. Abdesselam et al., Test of lepton flavor universality in ${B\to K^{\ast}\ell^{+}\ell^{-}}$ decays at Belle, arXiv:1904.02440.
  * (140) LHCb Collaboration, R. Aaij et al., Angular analysis and differential branching fraction of the decay $B^{0}_{s}\to\phi\mu^{+}\mu^{-}$, JHEP 09 (2015) 179, [arXiv:1506.08777].
  * (141) BaBar Collaboration, J. Lees et al., Measurement of the $B\to X_{s}l^{+}l^{-}$ branching fraction and search for direct CP violation from a sum of exclusive final states, Phys. Rev. Lett. 112 (2014) 211802, [arXiv:1312.5364].
  * (142) Belle Collaboration, S. Wehle et al., Lepton-Flavor-Dependent Angular Analysis of $B\to K^{\ast}\ell^{+}\ell^{-}$, Phys. Rev. Lett. 118 (2017), no. 11 111801, [arXiv:1612.05014].
  * (143) LHCb Collaboration, R. Aaij et al., Angular analysis of the $B^{0}\to K^{*0}e^{+}e^{-}$ decay in the low-q2 region, JHEP 04 (2015) 064, [arXiv:1501.03038].
  * (144) CMS Collaboration, A. M. Sirunyan et al., Angular analysis of the decay B${}^{+}\to$ K${}^{+}\mu^{+}\mu^{-}$ in proton-proton collisions at $\sqrt{s}=$ 8 TeV, Phys. Rev. D 98 (2018), no. 11 112011, [arXiv:1806.00636].
  * (145) LHCb Collaboration, R. Aaij et al., Differential branching fraction and angular analysis of $\Lambda^{0}_{b}\rightarrow\Lambda\mu^{+}\mu^{-}$ decays, JHEP 06 (2015) 115, [arXiv:1503.07138]. [Erratum: JHEP 09, 145 (2018)].
  * (146) LHCb Collaboration, R. Aaij et al., Angular moments of the decay $\Lambda_{b}^{0}\rightarrow\Lambda\mu^{+}\mu^{-}$ at low hadronic recoil, JHEP 09 (2018) 146, [arXiv:1808.00264].
  * (147) LHCb Collaboration, R. Aaij et al., Test of lepton universality with $B^{0}\rightarrow K^{*0}\ell^{+}\ell^{-}$ decays, JHEP 08 (2017) 055, [arXiv:1705.05802].
  * (148) BaBar Collaboration, J. P. Lees et al., Measurement of Branching Fractions and Rate Asymmetries in the Rare Decays $B\to K^{(*)}l^{+}l^{-}$, Phys. Rev. D 86 (2012) 032012, [arXiv:1204.3933].


Plugging in the upper bound on on $n$ from the lemma statement to (30) gives
us
$|\mathbb{P}_{\alpha,\beta}(A_{n})-\mathbb{P}_{\alpha,\beta+\epsilon}(A_{n})|<\frac{1}{2}$.

However, suppose _neither_ inequality in the lemma statement holds, then by
the definition of $A_{n}$, we have that
$|\mathbb{P}_{\alpha,\beta}(A_{n})-\mathbb{P}_{\alpha,\beta+\epsilon}(A_{n})|\geq|\frac{3}{4}-\frac{1}{4}|=\frac{1}{2}$
for $n^{3/2}\leq\frac{\alpha m}{\epsilon}$, which is a contradiction. ∎

Above lemma holds for any two prices $\pi^{*}_{1}\neq\pi^{*}_{2}$. Readers
should think of $\pi^{*}_{1},\pi^{*}_{2}$ as the optimal prices for customer
$n+1$ in the $(\alpha,\beta)$ and $(\alpha,\beta+\epsilon)$ market
respectively. To reduce clutter in the following Corollary statement, let
$\pi^{*}_{1}=\pi^{*}(x,\alpha,\beta,T^{\prime}),\pi^{*}_{2}=\pi^{*}(x,\alpha,\beta+\varepsilon,T^{\prime})$,
and let $\pi^{*}_{M}(x)$ be the optimal price for market $M$ i.e.,
$\pi^{*}_{M}=\pi^{*}_{1}$ if $M=(\alpha,\beta)$ and $\pi^{*}_{M}=\pi^{*}_{2}$
otherwise.

######  Corollary E.1.

Consider market $M$ sampled uniformly at random from
$\\{(\alpha,\beta),(\alpha,\beta+\epsilon)\\}$, let $M^{\prime}$ be the other
market. Suppose $n$ is such that $n^{3/2}\leq\frac{\alpha m}{\epsilon}$. Let
$[(p_{1},I_{1}),\ldots,(p_{n},I_{n})]$ be a sequence of $n$ observations
generated from the market $M$, where $p_{i}\in\mathcal{F}_{i-1}$. Fix any
pricing algorithm that outputs $\pi$ based on the first $n$ observations. Then
for any $x\in[0,1)$, any $T^{\prime}$ such that $\pi^{*}_{1}\neq\pi^{*}_{2}$:

$\mathbb{P}\left(\left|\pi-\pi^{*}_{M^{\prime}}(x)\right|\leq\left|\pi-\pi^{*}_{M}(x)\right|\right)\geq\frac{1}{8}$

where the probability is taken both with respect to the randomness from the
stochastic arrival times, as well as the uniform sampling of the market
parameters.

###### Proof.

This directly follows from Lemma E.3. The extra factor of $1/2$ in the
probability comes from the fact that we randomly picked one market out of two.
∎

###  E.3 Step 3: Lipschitz bound on the optimal pricing policy

######  Lemma E.4.

For any remaining time $T\geq\frac{(1+\sqrt{2})e}{\alpha+\beta}$ and
$x\leq\frac{\alpha^{2}e}{4(\alpha+\beta)^{3}T}$

$\frac{\partial\pi^{*}(x,\alpha,\beta,T)}{\partial\beta}\leq\frac{-\alpha
e}{4(\alpha+\beta)^{3}T}$

###### Proof.

Differentiating both sides of (17) with respect to $\beta$ gives us

$\frac{1}{1-X^{*}_{T}(x)}\frac{\partial
X^{*}_{T}(x)}{\partial\beta}=\frac{TX^{*}_{T}(x)}{2\beta
TX^{*}_{T}(x)+(\alpha-\beta)T+e}$

We omit the initial state argument $x$ and denote $X^{*}_{T}=X^{*}_{T}(x)$ in
the following proof.

$\displaystyle\frac{\partial\pi^{*}(x,\alpha,\beta,T)}{\partial\beta}$
$\displaystyle=\frac{1}{1-X^{*}_{T}}\frac{\partial
X^{*}_{T}}{\partial\beta}-\frac{X^{*}_{T}}{\alpha+\beta
X^{*}_{T}}-\frac{\beta\frac{\partial X^{*}_{T}}{\partial\beta}}{\alpha+\beta
X^{*}_{T}}+\frac{x}{\alpha+\beta x}$ $\displaystyle=\frac{TX^{*}_{T}}{2\beta
TX^{*}_{T}+(\alpha-\beta)T+e}-\frac{X^{*}_{T}}{\alpha+\beta
X^{*}_{T}}-\frac{\beta}{\alpha+\beta
X^{*}_{T}}\frac{(1-X^{*}_{T})TX^{*}_{T}}{2\beta
TX^{*}_{T}+(\alpha-\beta)T+e}+\frac{x}{\alpha+\beta x}$
$\displaystyle=\frac{TX^{*}_{T}(\alpha+\beta X^{*}_{T})-X^{*}_{T}(2\beta
TX^{*}_{T}+(\alpha-\beta)T+e)-(1-X^{*}_{T})\beta TX^{*}_{T}}{(\alpha+\beta
X^{*}_{T})(2\beta TX^{*}_{T}+(\alpha-\beta)T+e)}+\frac{x}{\alpha+\beta x}$
$\displaystyle=\frac{-eX^{*}_{T}}{(\alpha+\beta X^{*}_{T})(2\beta
TX^{*}_{T}+(\alpha-\beta)T+e)}+\frac{x}{\alpha+\beta x}$
$\displaystyle=\frac{-eX^{*}_{T}}{(\alpha+\beta
X^{*}_{T})\sqrt{(\alpha+\beta)^{2}T^{2}+2e(\alpha-\beta)T+e^{2}+4e\beta
xT}}+\frac{x}{\alpha+\beta x}$
$\displaystyle\leq\frac{-eX^{*}_{T}}{(\alpha+\beta
X^{*}_{T})((\alpha+\beta)T+e)}+\frac{x}{\alpha+\beta x}$

We replaced $X^{*}_{T}$ with (18) in the last equality. The last inequality
follows from the fact that
$\sqrt{(\alpha+\beta)^{2}T^{2}+2e(\alpha-\beta)T+e^{2}+4e\beta
xT}\leq(\alpha+\beta)T+e$. In particular, if
$T\geq\frac{(1+\sqrt{2})e}{\alpha+\beta}$, then
$\sqrt{(\alpha+\beta)^{2}T^{2}+2e(\alpha-\beta)T+e^{2}+4e\beta
xT}\leq\sqrt{2}(\alpha+\beta)T$. Also

$X^{*}_{T}(x,\alpha,\beta)\geq X^{*}_{T}(0,\alpha,\beta)\geq
X^{*}_{T}(0,\alpha,0)=\frac{\alpha T}{\alpha
T+e}\geq\frac{1+\sqrt{2}}{2+\sqrt{2}}\frac{\alpha}{\alpha+\beta}$

The first inequality is easy to see from (18), the second follows from
Corollary B.1, and the last inequality follows from the assumption on $T$. So
the above can be simplified to

$\frac{\partial\pi^{*}(x,\alpha,\beta,T)}{\partial\beta}\leq\frac{-e\frac{1+\sqrt{2}}{2+\sqrt{2}}\frac{\alpha}{\alpha+\beta}}{(\alpha+\beta)^{2}T\sqrt{2}}+\frac{x}{\alpha+\beta
x}=\frac{-\alpha e}{2(\alpha+\beta)^{3}T}+\frac{x}{\alpha+\beta x}$

Then, it easy to verify that for
$x\leq\frac{\alpha^{2}e}{4(\alpha+\beta)^{3}T}$,
$\frac{\partial\pi^{*}(x,\alpha,\beta,T)}{\partial\beta}\leq\frac{-\alpha
e}{4(\alpha+\beta)^{3}T}$ ∎

###  E.4 Step 4: Putting it all together for proof of Theorem 5.1

We are now ready to prove our main lower bound result Theorem 5.1. In the
following proof, as defined earlier in Step 1, $p_{x}$ denotes the induced
price trajectory from the algorithm’s offered prices, and $t^{det}_{x}$ is the
time when adoption level hits $x$ in the deterministic Bass model on following
$p_{x}$ (both were defined in detail in the paragraphs before Lemma E.1).

###### Proof of Theorem 5.1.

Let $T=\frac{2(1+\sqrt{2})}{\alpha+\beta}$,
$\epsilon=\frac{(\alpha+\beta)^{2}}{\alpha}$,
$n=\left\lfloor\left(\frac{\alpha}{\alpha+\beta}\right)^{4/3}m^{2/3}\right\rfloor$.
Randomly draw the Bass model parameters from
$\\{(\alpha,\beta),(\alpha,\beta+\epsilon)\\}$ with equal probabilities. We
denote the chosen market as $M$, and the other market $M^{\prime}$. In the
calculations that follow, we use $\mathbb{E}_{M}$ to indicate that the
expectation is taken with respect to both the random choice of Bass model
parameters as well as the randomness in the stochastic Bass model.

To reduce clutter, let
$\pi^{*}_{1}(x)=\pi^{*}(x,\alpha,\beta,T-t^{det}_{x}),\pi^{*}_{2}(x)=\pi^{*}(x,\alpha,\beta+\varepsilon,T-t^{det}_{x})$,
and let $\pi^{*}_{M}(x)$ be the optimal price i.e.,
$\pi^{*}_{M}(x)=\pi^{*}_{1}(x)$ if $M=(\alpha,\beta)$ and
$\pi^{*}_{M}(x)=\pi^{*}_{2}(x)$ otherwise.

$\displaystyle\mathbb{E}_{M}[\text{Pseudo-Regret}(T)]$
$\displaystyle\text{(Lemma~{}\ref{lem: lower bound pseudo regret price
difference})}\geq$
$\displaystyle\mathbb{E}_{M}\left[m\int_{0}^{n/m}\min\left(\frac{1}{4}(\pi^{*}_{M}(x)-p_{x})^{2},\frac{1}{10}\right)dx\right]-O\left(\frac{(\alpha+\beta)^{1/3}Te^{C}}{\alpha^{1/3}}m^{1/3}\sqrt{\log(m)}\right)$
$\displaystyle\geq$ $\displaystyle
m\int_{0}^{n/m}\mathbb{E}_{M}\left[\min\left(\frac{1}{4}\left(\frac{\pi^{*}_{1}(x)-\pi^{*}_{2}(x)}{2}\right)^{2},\frac{1}{10}\right)\mathbbm{1}\left(\left|p_{x}-\pi^{*}_{M^{\prime}}(x)\right|\leq\left|p_{x}-\pi^{*}_{M}(x)\right|\right)\right]dx$
$\displaystyle-O\left(\frac{m^{1/3}e^{C}}{(\alpha+\beta)^{2/3}\alpha^{1/3}}\sqrt{\log(m)}\right)$
$\displaystyle\text{(Lemma~{}\ref{lem: optimal price difference
deterministic})}\geq$ $\displaystyle
m\int_{0}^{n/m}\mathbb{E}_{M}\left[\min\left(\left(\frac{\alpha
e\varepsilon}{16(\alpha+\beta)^{3}(T-t^{det}_{x})}\right)^{2},\frac{1}{10}\right)\mathbbm{1}\left(\left|p_{x}-\pi^{*}_{M^{\prime}}(x)\right|\leq\left|p_{x}-\pi^{*}_{M}(x)\right|\right)\right]dx$
$\displaystyle-O\left(\frac{m^{1/3}e^{C}}{(\alpha+\beta)^{2/3}\alpha^{1/3}}\sqrt{\log(m)}\right)$
$\displaystyle\text{($T\geq T-t^{det}_{x}$)}\geq$ $\displaystyle
m\int_{0}^{n/m}\min\left(\left(\frac{e}{32(1+\sqrt{2})}\right)^{2},\frac{1}{10}\right)\mathbb{E}_{M}\left[\mathbbm{1}\left(\left|p_{x}-\pi^{*}_{M^{\prime}}(x)\right|\leq\left|p_{x}-\pi^{*}_{M}(x)\right|\right)\right]dx$
$\displaystyle-O\left(\frac{m^{1/3}e^{C}}{(\alpha+\beta)^{2/3}\alpha^{1/3}}\sqrt{\log(m)}\right)$
$\displaystyle\text{(Corollary~{}\ref{cor:decision rule limit})}\geq$
$\displaystyle\frac{n}{8}\min\left(\left(\frac{e}{32(1+\sqrt{2})}\right)^{2},\frac{1}{10}\right)-O\left(\frac{m^{1/3}e^{C}}{(\alpha+\beta)^{2/3}\alpha^{1/3}}\sqrt{\log(m)}\right)$
$\displaystyle=$
$\displaystyle\Omega\left(\left(\frac{\alpha}{\alpha+\beta}\right)^{4/3}m^{2/3}\right)-O\left(\frac{m^{1/3}e^{C}}{(\alpha+\beta)^{2/3}\alpha^{1/3}}\sqrt{\log(m)}\right)$
$\displaystyle=$
$\displaystyle\Omega\left(\left(\frac{\alpha}{\alpha+\beta}\right)^{4/3}m^{2/3}\right)$

Using the fact that the maximum expected regret between the two markets must
be at least the average, we have that for at least one of $(\alpha,\beta)$ and
$(\alpha,\beta+\epsilon)$,

$\displaystyle\mathbb{E}[\text{Pseudo-
Regret}(T)]\geq\Omega\left(\left(\frac{\alpha}{\alpha+\beta}\right)^{4/3}m^{2/3}\right)$

Finally, we still need to check the assumptions of Lemma E.2, Lemma E.4 and
Corollary E.1 are satisfied for large enough $m$.

To apply Lemma E.2, we need
$X^{*}_{T}\geq\left(\frac{\alpha}{\alpha+\beta}\right)^{4/3}m^{-1/3}=\frac{n}{m}$.
Using the expanded notation $X^{*}(x,\alpha,\beta)$ introduced in Appendix B,
we know from Corollary B.1 that $X^{*}_{T}(0,\alpha,\beta)\geq
X^{*}_{T}(0,\alpha,0)=\frac{\alpha T}{\alpha T+e}$. It is easy to verify that
for large enough $m$, specifically when
$m^{1/3}\geq\left(\frac{\alpha}{\alpha+\beta}\right)^{4/3}\left(2+\frac{e(\alpha+\beta)}{\alpha(1+\sqrt{2})}\right)$
and $T=2(1+\sqrt{2})e/(\alpha+\beta)$, $\frac{n}{m}\leq\frac{\alpha T}{\alpha
T+e}\leq X^{*}_{T}(0,\alpha,\beta)=X^{*}_{T}$.

To apply Lemma E.4 we need the conditions $T^{\prime}_{x}\coloneqq
T-t^{det}_{x}\geq\frac{(1+\sqrt{2})e}{\alpha+\beta}$ and
$x\leq\frac{\alpha^{2}e}{4(\alpha+\beta)^{3}T^{\prime}_{x}}$ for all
$x\leq\frac{n}{m}$. For any fixed $\alpha,\beta$, for large enough $m$,
specifically for for
$m^{1/3}\geq\frac{\alpha^{1/3}}{(\alpha+\beta)^{4/3}}(\alpha+e^{C})$, where
$C$ is a function of $\alpha,\beta$ as defined in Assumption 2, we have that
$t^{det}_{n/m}\leq\frac{n}{e^{-p_{max}}\alpha(m-n)}=\frac{e^{p_{max}}n}{\alpha(m-n)}\leq\frac{(1+\sqrt{2})e}{\alpha+\beta}$.
The last inequality followed from simple algebraic calculations using the
assumption on $m$ and Assumption 2. Therefore,
$T-t^{det}_{x}\geq\frac{(1+\sqrt{2})e}{\alpha+\beta}$ for all
$x\leq\frac{n}{m}$. This satisfies the first condition. Furthermore, for
$m^{1/3}\geq 8(1+\sqrt{2})\left(\frac{\alpha+\beta}{\alpha}\right)^{2/3}$,
$x\leq
n/m=\left(\frac{\alpha}{\alpha+\beta}\right)^{4/3}m^{-1/3}\leq\frac{\alpha^{2}}{8(1+\sqrt{2})(\alpha+\beta)^{2}}=\frac{\alpha^{2}e}{4(\alpha+\beta)^{3}T}\leq\frac{\alpha^{2}e}{4(\alpha+\beta)^{3}T^{\prime}_{x}}$.
This satisfies the second condition.

Finally, to apply Corollary E.1 we needed the condition
$n^{3/2}\leq\frac{\alpha m}{\epsilon}$. By plugging in our choice of
$m,\epsilon$, it is easy to verify that this condition is satisfied.

∎

##  Appendix F Auxiliary Results

######  Lemma F.1.

Let $\bar{p}$ be a constant upper bound on the prices that the seller can
offer: $p_{t}\leq\bar{p}$. Then there exists constants $C_{0},C_{1}$
independent of $m$ such that for $T\geq C_{0}\log(m)+C_{1}$, there is a
trivial policy that achieves $O\left(\log(\frac{1}{\delta})\right)$ regret.

###### Proof.

Consider the trivial pricing policy where the seller sets the price to the
highest possible value $\bar{p}$ for the entire planning horizon. The rate of
arrival of customer $d$ is given by
$\lambda(\bar{p},\frac{d-1}{m})=e^{-\bar{p}}(\alpha+\beta\frac{d-1}{m})(m-d+1)\geq
e^{-\bar{p}}\alpha(m-d+1)$

We now divide the market of $m$ customers into a sequence of segments with
geometrically decreasing lengths. Let $m_{0}=0,m_{1}=\lfloor
m/2\rfloor,m_{2}=m_{1}+\lfloor m/4\rfloor,\ldots,m_{i}=m_{i-1}+\lfloor
m/2^{i}\rfloor$. Let
$K=\lfloor\log_{2}(m)-\log_{2}\log(\frac{1}{\delta})-2\rfloor$ such that
$m_{K}=m-O(\log(\frac{1}{\delta}))$. We call customers $d$ for $m_{i-1}<d\leq
m_{i}$ the customers in segment $i$.

Since the price is constant, we use the following short hand notation
$\lambda_{d}\coloneqq\lambda(\bar{p},\frac{d-1}{m})$. Note that by
construction, $m-m_{i}\geq m/2^{i}$. So for $m_{i-1}<d\leq m_{i}$,
$\lambda_{d}\geq\underline{\lambda}_{i}\coloneqq
e^{-\bar{p}}\alpha\frac{m}{2^{i}}$.

Let $I_{d}$ be the stochastic inter-arrival time between customer $d-1$ and
$d$, then using Lemma A.1, we can obtain the following bound on the time it
takes for all customers in segment $i$ to arrive:

$\displaystyle\mathbb{P}\left(\sum\limits_{d=m_{i-1}+1}^{m_{i}}I_{d}-\frac{1}{\lambda_{d}}\geq\epsilon_{i}\right)\leq
exp\left(\frac{-\epsilon_{i}^{2}\underline{\lambda}_{i}^{2}}{8(m_{i}-m_{i-1})}\right)$

Setting this to $\delta$ and solving for $\epsilon_{i}$, one can verify that
with probability $1-\delta$,

$\sum\limits_{d=m_{i-1}+1}^{m_{i}}I_{d}-\frac{1}{\lambda_{d}}\geq\epsilon_{i}\coloneqq\frac{4e^{\bar{p}}}{\alpha}\sqrt{\log(\frac{1}{\delta})\frac{2^{i}}{m}}$

Note that we needed
$\epsilon_{i}\leq\frac{2(m_{i}-m_{i-1})}{\underline{\lambda}_{i}}$ to apply
Lemma A.1. This condition is satisfied for $i=1,\ldots,K$ because:

$\displaystyle\epsilon_{i}\leq\frac{2(m_{i}-m_{i-1})}{\underline{\lambda}_{i}}\iff$
$\displaystyle\sqrt{\log(\frac{1}{\delta})\frac{2^{i}}{m}}\leq\frac{1}{2}\iff
i\leq
2\log_{2}\left(\frac{1}{2}\sqrt{\frac{m}{\log(1/\delta)}}\right)=\log(m)-\log_{2}\log(\frac{1}{\delta})-2$

Applying a union bound to this result for all segments, we have that with
probability at least $1-\delta\log_{2}(m)$, we have that

$\displaystyle\sum\limits_{d=1}^{m_{K}}I_{d}\leq$
$\displaystyle\sum\limits_{d=1}^{m_{K}}\frac{1}{\lambda_{d}}+\sum_{i=1}^{K}\epsilon_{i}$
$\displaystyle\leq$
$\displaystyle\sum_{i=1}^{K}\frac{m_{i}-m_{i-1}}{\underline{\lambda}_{i}}+\sum_{i=1}^{K}\frac{4e^{\bar{p}}}{\alpha}\sqrt{\log(\frac{1}{\delta})\frac{2^{i}}{m}}$
$\displaystyle=$ $\displaystyle
2K\frac{e^{\bar{p}}}{\alpha}+\frac{4e^{\bar{p}}}{\alpha}\sqrt{\log(\frac{1}{\delta})\frac{1}{m}}\frac{\sqrt{2}(2^{K/2}-1)}{\sqrt{2}-1}$
$\displaystyle\leq$
$\displaystyle\frac{e^{\bar{p}}}{\alpha}\left(2\log_{2}(m)+\frac{4\sqrt{2}}{\sqrt{2}-1}\sqrt{\log(\frac{1}{\delta})}\right)$
$\displaystyle=$ $\displaystyle C_{0}\log(m)+C_{1}$

This means that if $T\geq C_{0}\log(m)+C_{1}$, then with probability
$1-\delta\log_{2}(m)$, the realized revenue is at least $m_{K}\bar{p}$. Since
the maximum possible revenue is $m\bar{p}$, the regret is at most
$\bar{p}(m-m_{K})=O(\log(\frac{1}{\delta}))$. ∎