Pierre-O Goffard

Sequential Monte Carlo samplers to fit and compare insurance loss models

Keywords: 60G55, 60G40, 12E10 Composite model, Bayesian statistics, Sequential Monte Carlo sampler

Insurance loss distributions are characterized by a high frequency of small claim amounts and a lower, but not insignificant, occurrence of large claim amounts. Composite models, which link two probability distributions, one for the "body" and the other for the "tail" of the loss distribution, have emerged in the actuarial literature to take this specificity into account. The parameters of these models summarize the distribution of the losses. One of them corresponds to the breaking point between small and large claim amounts. The composite models are usually fitted using maximum likelihood estimation. A Bayesian approach is considered in this work. Sequential Monte Carlo samplers are used to sample from the posterior distribution and compute the posterior model evidences to both fit and compare the competing models. The method is validated via a simulation study and illustrated on an insurance loss dataset.

Introduction

The distribution of losses in property and casualty insurance is characterized by a high frequency of small claim amounts and a lower, but not insignificant, frequency of considerably larger claim amounts. Composite models, that combine two models one for the body and the other for the tail of the loss distribution, have emerged in the actuarial science litterature as a response to this specificity. A model that accounts for both 1 extreme and moderate claim amounts is necessary to inform the insurance company decision process regarding premium calculation, risk capital allocation and risk transfer optimization.

The breakpoint between small and large claims is a key parameter in composite models. Two approaches are commonly used in practice to determine it. The first, known as the fixed-threshold approach, consists of choosing a threshold using tools from extreme value theory before adjusting the components of the composite models. It relies on famous graphical visualizations such as the mean-excess plot [START_REF] Ghosh | A discussion on mean excess plots[END_REF], the Hill plot [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF] or the Gerstengarbe plot [START_REF] Gerstengarbe | A method for the statistical definition of extreme-value regions and their application to meteorological time series[END_REF]. Automatic threshold selection methods have been implemented to mitigate the subjectivity inherent in locating a stable area on a curve. The goal is often to minimize the root mean squared error around the tail index estimate, see Caeiro and Gomes [START_REF] Caeiro | Extreme value modeling and risk analysis : methods and applications, chapter Threshold selection in extreme value analysis[END_REF], or to pass a goodness-of-fit test, see Guillou and Hall [START_REF] Guillou | A diagnostic for selecting the threshold in extreme value analysis[END_REF]. These techniques rely on the hypothesis that the data points in excess of the threshold are distributed according to a Generalized Pareto Distribution (GPD). For an overview on extreme value data analysis, I refer the reader to the textbook of Beirlant et al. [START_REF] Beirlant | Statistics of Extremes[END_REF]. Another approach, less used in practice but widely studied in the actuarial literature, treats the threshold like any other parameter of the composite model and performs a simultaneous estimation. I will refer to this procedure as the free-threshold approach. A survey on the use of composite models and threshold selection methods on insurance data is conducted in the work of Wang et al. [START_REF] Wang | Modelling extreme claims via composite models and threshold selection methods[END_REF]. Maximum likelihood estimation is commonly used to fit several combinations of models for the body and the tail of the claim amounts distribution. The adequacy of each model is then measured using standard information criteria, see for instance Grün and Miljkovic [START_REF] Grün | Extending composite loss models using a general framework of advanced computational tools[END_REF].

The present work proposes to fit and compare composite models in a Bayesian way. Bayesian statistics take the model parameters to be random variables. Inference is drawn from the posterior distribution of the parameters obtained by updating the a priori assumptions via the likelihood function of the data, for an overview see the book by Geldman et al. [START_REF] Gelman | Bayesian Data Analysis[END_REF]. The posterior distribution is often unavailable and must be approximated by an empirical distribution. Markov chain Monte Carlo (MCMC) simulation schemes, such as the well known Metropolis-Hasting and Gibbs samplers, have become the go to techniques to sample from the posterior distribution. Probabilistic programming softwares like WINBUGS [START_REF] Lunn | Winbugs-a bayesian modelling framework: concepts, structure, and extensibility[END_REF], JAGS [START_REF] Plummer | A program for analysis of bayesian graphical models using gibbs sampling[END_REF], STAN [START_REF] Carpenter | Stan: A probabilistic programming language[END_REF] and PYMC [START_REF] Salvatier | Probabilistic programming in python using PyMC3[END_REF], have been designed over the years so that practitioners do not have to worry about the fine tuning of these sophisticated algorithms. The lognormal-Pareto composite model of Scollnik [START_REF] Scollnik | On composite lognormal-pareto models[END_REF] is actually an example of the WINBUGS documentation, see [START_REF] Lunn | Winbugs-a bayesian modelling framework: concepts, structure, and extensibility[END_REF]Examples Vol. III]. The application of Bayesian techniques to composite models enables to quantify the uncertainty around the threshold parameter, as noted in the survey of Scarrot and MacDonald [START_REF] Scarrott | A review of extreme value threshold es-timation and uncertainty quantification[END_REF]. Instead of the standard MCMC sampler, a Sequential Monte Carlo Sampler (SMC) is put together. This algorithm builds a sequence of empirical distributions, made of weighted particles, that targets the posterior distribution during the final iteration. Generic SMC samplers are described in the seminal paper of Del Moral et al. [START_REF] Del Moral | Sequential monte carlo samplers[END_REF]. An approximation of the posterior distribution normalizing constant, referred to as the marginal likelihood, follows from the weights of the successive particle clouds. MCMC algorithms bypass the evaluation of this constant which is nevertheless necessary for the evaluation of Bayes factors to select the right model, see Kass and Raftery [START_REF] Kass | Bayes factors[END_REF]. In addition to providing an approximation of the marginal likelihood, SMC samplers can sample from complicated multimodal posterior distributions, save the trouble of tuning some hyperparameters and are easy to paralellize which is a key feature in the era of multi-core processor computers. The multimodality of the posterior distribution around the threshold parameter has been encountered in practice by Cabras and Castellanos [8, Figure 3].

The main contribution of this work is to provide an efficient computational tool to fit and compare a myriad of composite models with a fixed or free threshold. This tool takes the form of a Python package available on pip1 for anyone to use. The posterior distribution of the parameters and the posterior probability of the models make it possible to account for the uncertainty around the values of the parameters and the model to be used. This uncertainty can then be taken into account when estimating quantities of interest to risk managers, such as extreme quantiles or probabilities of ruin at one year, via credible sets. It is also possible to use the posterior probabilities of the model as weights to combine the estimate provided by each model. This procedure, known as Bayesian Model Averaging (BMA), see Kass and Raftery [START_REF] Kass | Bayes factors[END_REF], is investigated. This work can be viewed as the Bayesian counterpart of the work of Grün and Miljkovic [START_REF] Grün | Extending composite loss models using a general framework of advanced computational tools[END_REF]. Besides the inference method, the main difference lies in the regularity assumptions of the probability density function (pdf) at the threshold parameter. Namely, the authors of [START_REF] Grün | Extending composite loss models using a general framework of advanced computational tools[END_REF] impose continuity and differentiability at the threshold whereas the models considered here can be discontinuous or only continuous but not differentiable. Another difference is that the algorithm allows the user to set in advance or not the breaking point between low and high severities to enable the comparison between the fixed and free threshold approach.

The remainder of the paper is organized as follows. Section 2 provides a brief overview on composite loss models and their use for risk management purposes. Section 3 recalls the principles of Bayesian statistics and presents the algorithmic details of the sequential Monte Carlo samplers used to fit the composite models. A simulation experiment is conducted in Section 4 to assess the consistency and finite-sample performance of the estimation and model selection procedures. Section 5 illustrates the application of the SMC algorithm on Losses in insurance are usually modelled by nonnegative random variables with probability density function (pdf) denoted by f (x; θ), where θ ∈ Θ ⊂ R d is the parameter space and d > 0 its dimension. The goal is then to find the parameter value θ that best explains the data x = (x 1 , . . . , x n ). This is usually achieved by maximizing the likelihood function L(x|θ). In the case of independent and identically distributed (iid) data (which is the case considered here), the likelihood function is given by

L(x|θ) = n i=1 f (x i ; θ).
The occurence of extreme losses makes the vast majority of the simple parametric models (2-3 parameters) inadequate. The lack of fit of these simple models leads to consider more flexible ones.

Definition and assumptions

Composite models, sometimes called splicing or spliced models, combine two models, one for the "body" and the other for the "tail" of the distribution. The pdf of a composite model is defined as

f (x) =            p f 1 (x) F 1 (γ) , si x < γ, (1 -p) f 2 (x) 1-F 2 (γ) , si x ≥ γ, (1)
where f 1 , F 1 , f 2 , and F 2 are the pdf and cumulative distribution function (cdf) of the body and the tail of the loss distribution respectively. The parameter p ∈ (0, 1) is referred to as the mixing parameter. The threshold parameter γ > 0 is the breaking point that distinguishes small claim amounts from large ones. Such models appeared in the statistical literature as early as the work of Mendes and Lopez [START_REF] Vaz De Melo Mendes | Data driven estimates for mixtures[END_REF] who used maximum likelihood estimation to infer the parameters. A Bayesian framework was proposed by Berhens et al. [START_REF] Behrens | Bayesian analysis of extreme events with threshold estimation[END_REF] to enable the quantification of the uncertainty around the threshold parameter. A commonly used simplification of model [START_REF] Abu Bakar | Modeling loss data using composite models[END_REF]

sets p = F 1 (γ), (2)
making γ the p-quantile of the F 1 distribution. It has been pointed out several times that the potential discontinuity at x = γ of the density (1) could pose problems in some practical situations. Continuity at the

threshold f (γ -) = f (γ + ) is achieved by setting p = f 1 (γ) F 1 (γ) f 2 (γ) 1 -F 2 (γ) . ( 3
)
authors further imposed differentiability at the threshold f ′ (γ -) = f ′ (γ + ) which leads to loose another degree of freedom by letting γ be the solution of

f ′ 1 (γ) f 1 (γ) - f ′ 2 (γ) f 2 (γ) = 0, (4)
see the work of Grün and Miljkovic [START_REF] Grün | Extending composite loss models using a general framework of advanced computational tools[END_REF]. Wang et al [START_REF] Wang | Modelling extreme claims via composite models and threshold selection methods[END_REF] have reported that imposing such a condition tends to make the model not flexible enough leading to the conclusion that practitionners would be better off by choosing a threshold before fitting a general composite model.

In the present paper, the bulk and tail of composite models are selected from the list in Table 1. In contrast to Name Parameters pdf Remark 2.1. Imposing differentiability at x = γ entails the loss of two degrees of freedom by letting the threshold be the solution of (4) and therefore depend on the other parameters of the composite models. The practical issue is that the uniqueness of the solutions of (4) is far from being granted. A numerical root finding procedure would be necessary which will slow down the inference process. A tedious case by case study of the 120 possible composite models would be necessary. Such a study has been carried out in the litterature for a few particular cases including the lognormal-Pareto case, see Cooray and Ananda [START_REF] Cooray | Modeling actuarial data with a composite lognormal-pareto model[END_REF] and Scollnik [START_REF] Scollnik | On composite lognormal-pareto models[END_REF], and the Weibull-Pareto case, see Scollnik and Sun [START_REF] Scollnik | Modeling with weibull-pareto models[END_REF] and Abu Bakar et al. [START_REF] Abu Bakar | Modeling loss data using composite models[END_REF]. In addition to general composite model being discontinuous or continuous, The

Exponential Exp(δ) δ > 0 δe -δx , x > 0 Gamma Gamma(r, m) r, m > 0 x r-1 e -x/m Γ (r)m r , x > 0 Weibull Weibull(k, β) k, β > 0 k β x β k-1 e -(x/β) k , x > 0 Lognormal Lognormal(µ, σ ) µ ∈ R, σ > 0 1 xσ √ 2π exp - (log(x)-µ) 2 2σ 2 , x > 0 Inverse-Gaussian Inverse-Gaussian(µ, λ) µ, λ > 0 λ 2πx 3 exp λ(x-µ) 2 ) 2µ 2 x , x > 0 Inverse-Gamma Inverse-Gamma(r, m) r, m > 0 e -m/x m r x r+1 γ(r) , x > 0 Inverse-Weibull Inverse-Weibull(k, β) k, β > 0 kβ k x -k-1 e -(β/x) k , x > 0 Lomax Lomax(α, σ ) α, σ > 0 ασ α (σ +x) α+1 , x > 0 Log-Logistic Log-Logistic(β, σ ) β, σ > 0 βσ β x β-1 (σ β +x β ) 2 , x > 0 Burr Burr(α, β, σ ) α, β, σ > 0 αβσ αβ x β-1 (σ β +x β ) α+1 , x > 0 Pareto Pareto(α, γ) α, γ > 0 γ α x α+1 , x > γ Generalized Pareto GPD(ξ, σ , γ) ξ, σ , γ > 0 σ -1 1 + ξ(x-γ) σ -(ξ+1)/ξ , x ≥ 0
m = γ k + α , p = αΓ (k)F 1 (γ; r, m)e k+α (k + α) -k 1 + αΓ (k)F 1 (γ; r, m)e k+α (k + α) -k ,
where F 1 (γ; r, m) is the gamma distribution cdf.

The bayes-splicing package is able to fit a total of 375 loss models when accounting for all the distributions in Table 1 either considered separately or as component of a splicing model (with three possible settings "discontinuous", "simple" or "continuous"). The models are fitted and compared using state of the art Bayesian statistics computational tools which are presented in Section 3. The goal is to provide the right take on the loss distribution to assess the level of capital required for solvency purposes defined in Section 2.2.

Composite models usage in actuarial science

Over a given time period, a year say, a non-life insurance company handles a random number of claims N , each of which is associated to a randomly sized compensation forming hereby a sequence U 1 , . . . , U N . Usually, the U ′ i s are taken to be independent and identically distributed random variables, independent from the claim frequency N . The total cost over one year amounts to

S = N i=1 U i .
(

) 5
The random sum (5) corresponds to the liability of the insurer within what is called a collective model, see for instance the book of Klugman et al. [START_REF] Klugman | Loss Models: From Data to Decisions[END_REF]. Actuaries and risk managers typically want to quantify the risk of large losses by a single comprehensible number, a risk measure. One popular risk measure is the Value-at-Risk which corresponds to a high order level quantile of the distribution of S. It is used by risk managers in banks, insurance companies, and other financial institutions to allocate risk reserves and to determine risk margins.

The right tail of U has a strong influence on the higher order quantile of S and must be modelled appropriately, using for instance a composite model. One way to mitigate the risk associated with large claims is to transfer some of it to a reinsurance company. The aggregate claim sizes is then divided into

S = D + R,
where D is the amounts that stays with the first-line insurer after reinsurance and R is the amount paid by the reinsurer. The first-line insurer receives premiums, that sum up to Π, from the policyholders and part of it is ceded to the reinsuer. We have

Π = Π R + Π D ,
where Π R is the reinsurance premium ceded to the reinsurer while Π D is the premium retained by the first-line insurer. A common reinsurance agreement in property and casualty insurance is the Excess-of-Loss (xol) one with

R = N i=1 min{(U i -P ) + , L}, and D = N i=1 min{U i , P }I {U i ≤P +L} + (U i -L)I {U i >P +L} , (6)
where (•) + denotes the positive part, P is the priority and L is the limit. The first-line insurer seeks values of P and L to optimize its expected surplus under some solvency constraints. We focus in this work on the left quantile of the distribution of Π D -D to determine a risk reserve. For a comprehensive overview on the statistical and actuarial aspect of reinsurance, the reader is referred to the book of Albrecher et al. [2].

The claim sizes are assumed to be iid from a composite model. The claim frequency is Poisson distributed with mean λ > 0. The premiums are defined by the expectation principle with a safety loading of 5% and computed using numerical integration2 . Once the parameters of the composite model have been inferred

then the quantiles of the insurer's surplus are calculated via crude Monte Carlo simulations. Through this application, we see the importance of precisely modeling the entire distribution of losses, which is the point of composite models. They are essential for an accurate assessment of insurance and reinsurance premiums but also of the probability of occurrence of serious claims. Let us move on to the Bayesian estimation procedure to fit and compare the composite models.

Bayesian inference via a sequential Monte Carlo sampler

Bayesian statistics defines the posterior distribution of the model parameters θ given the data x = (x 1 , . . . , x n ) as

π(θ|x) = L(x|θ)π(θ) Z(x) . ( 7
)
The posterior distribution [START_REF] Bladt | Threshold selection and trimming in extremes[END_REF] follows from applying Bayes' rule to update the prior distribution π(θ) using the likelihood function L(x|θ). Credible sets as well as point estimates can then be derived from π(θ|x) to draw inference on θ. The only issue is the denominator in [START_REF] Bladt | Threshold selection and trimming in extremes[END_REF] which is a normalizing constant given by

Z(x) = Θ L(x|θ)π(θ)dθ. ( 8
)
The above integral rarely admits a closed-form expression except when the model for the data has a conjugate prior distribution. Unfortunately, conjugate prior distributions almost only arise in exponential families of probability distributions, see Diaconis and Ylvisaker [START_REF] Diaconis | Conjugate priors for exponential families[END_REF]. In practice, one samples from the posterior distribution via Markov Chain Monte Carlo (mcmc) schemes. The Metropolis-Hasting random walk builds a sequence (θ i ) i≥0 by applying a Markov kernel K H (•|θ i ) to the current parameter value θ i , i ≥ 0. The parameter H corresponds to the magnitude of the perturbation. A new parameter value

θ * ∼ K H (•|θ i ) is accepted with probability α(θ i , θ * ) = max 1, L(x|θ * )π(θ * )K H (θ * |θ i ) L(x|θ 1 )π(θ i )K H (θ i |θ * ) , (9)
in which case θ i+1 = θ * , otherwise θ i+1 = θ i . The resulting sequence (θ i ) i≥0 forms a Markov chain trajectory having the posterior distribution as limiting distribution. A standard choice for the Markov kernel is the multivariate normal distribution

K H (•|θ) ∼ Norm(µ = θ, Σ = H), ( 10
)
where H is a matrix that matches the dimension of θ. The Metropolis-Hasting random walk efficiency, understood as the speed of convergence of the Markov chain toward its asymptotic distribution, decreases with the dimension of the parameters. The workaround is to turn to another well-known mcmc technique that generates a sequence (θ i ) i≥0 called Gibbs sampling. To sample from a multivariate posterior distribution π(θ|x), a Gibbs sampler samples from the univariate conditional distributions defined as

π(θ j |x, θ 1 , . . . , θ j-1 , θ j+1 , . . . , θ d ), for j = 1, . . . , d,
where θ = (θ 1 , . . . , θ d ). The current parameter value θ i is updated component per component starting with the first one

θ i 1 ∼ π(•|x, θ i 2 , . . . , θ i d ), (11)
before moving to the second one

θ i 2 ∼ π(•|x, θ i 1 , θ i 3 , . . . , θ i d ), (12)
and so on. The sequence (θ i ) i≥0 forms a Markov chain trajectory whose limiting distribution is the posterior distribution π(θ|x). The marginal distributions in [START_REF] Cooray | Modeling actuarial data with a composite lognormal-pareto model[END_REF], [START_REF] Danielsson | Tail index estimation: Quantile driven threshold selection[END_REF], etc., are usually unknown and one actually uses the Metropolis-Hasting scheme within the Gibbs iterations to sample from them. The Metropolis-Hasting within Gibbs type algorithms admit some drawbacks. First the algorithm must be initialized. In practice, several chains are launched from different starting points θ 0 to verify if they all converge toward the same distribution.

Second, the H parameter of the multivariate normal kernel in [START_REF] Carpenter | Stan: A probabilistic programming language[END_REF] must be tuned to ensure good sampling properties. It should reflect the variance of the posterior distribution which is unknown. The practical solution is to implement an adaptive procedure to adjust H on the fly to reach an acceptance rate of 23.4% which is deemed optimal, see the work of Roberts et al. [START_REF] Roberts | Weak convergence and optimal scaling of random walk metropolis algorithms[END_REF]. Third, the trajectory generation cannot be parallelized.

Lastly, mcmc algorithms allow one to sample the posterior distribution of any models as long as the likelihood function has a tractable expression by avoiding the evaluation of the normalizing constant in [START_REF] Bladt | Threshold selection and trimming in extremes[END_REF]. Indeed, the latter does not appear in the acceptance probability expression in [START_REF] Caeiro | Extreme value modeling and risk analysis : methods and applications, chapter Threshold selection in extreme value analysis[END_REF]. The normalizing constant is nevertheless important for Bayesian model selection as explained below.

Consider a set of competing models M = {m 1 , . . . , m J } and define a random variable M having a Probability Mass Function (pmf) concentrated on M. A prior distribution such that P(M = m j ) = π(m j ) ≥ 0, for j = 1, . . . , J, and J j=1 π(m j ) = 1 can then be specified and updated given the data to yield the posterior model evidence as

π(m j |x) = L(x|m j )π(m j ) J i=1 L(x|m i )π(m i ) , j = 1, . . . , J. ( 13
)
The likelihood L(x|m) of model m ∈ M follows from integrating over the possible values of the parameter θ as

L(x|m) = Θ L(x|m, θ)π(θ|m)dθ,
which corresponds exactly to the normalizing constant in [START_REF] Cabras | A bayesian approach for estimating extreme quantiles under a semiparametric mixture model[END_REF]. The best model achieves the highest model evidence [START_REF] Davison | Statistical Models[END_REF]. Another use of the posterior probabilities is the weighting of the estimation of the different models. Denote by ∆ a quantity. It is possible to combine the estimation ∆ j of each model m j as

∆ = J j=1 ∆ j π(m j |x).
This ensemble estimation procedure, known as Bayesian Model Averaging, is detailed in the work of Hoeting et al. [START_REF] Hoeting | Bayesian model averaging: A tutorial[END_REF] and tested out in Section 4.

This section describes a sequential Monte Carlo algorithm which allows one to sample from any posterior distributions while providing an approximation of the normalization constant. The implementation is effortless to parallelize and its hyperparameters are straigtforward to tune. Section 3.1 provides a quick reminder of the importance sampling principle required to understand the smc algorithm detailed in Section 3.2.

Importance sampling

Bayesian inference reduces to evaluating quantities such as

E π(θ|x) (ϕ) = Θ ϕ(θ)π(θ|x)dθ, ( 14
)
where E π(θ|x) is the expectation operator with respect to the posterior distribution and ϕ is some measurable application. The posterior mean, often used as point estimate, corresponds to the case ϕ(θ) = θ. The expectation ( 14) is evaluated through its Monte Carlo approximation

E π(θ|x) (ϕ) ≈ 1 K K i=1 ϕ(θ i ), (15)
where θ 1 , . . . , θ K is an iid sample distributed as π(θ|x). Importance sampling samples from a distribution g, instead of π(θ|x), either because it is more convenient or because that reduces the variance associated to the Monte Carlo estimator [START_REF] Diaconis | Conjugate priors for exponential families[END_REF]. The approximation of the normalizing constant relies on the following identity

E π(θ|x) (ϕ) = Θ ϕ(θ)π(θ|x)dθ = Θ ϕ(θ) L(x|θ)π(θ) Z(x) dθ = Z(x) -1 Θ ϕ(θ) L(x|θ)π(θ) g(θ) g(θ)dθ = Z(x) -1 Θ ϕ(θ)w(θ)g(θ)dθ = Z(x) -1 E g (ϕ • w),
where w(θ) = L(x|θ)π(θ)/g(θ) is an unnormalized weight function. Taking ϕ(θ) = 1 yields the following expression of the normalizing constant

Z(x) = E g (w),
which may be approximated by

Z(x) ≈ 1 K K i=1 w( θ i ),
where θ 1 , . . . , θ N is an iid sample generated from the proposal g. Importance sampling ultimately yields a cloud of weighted particles {W i , θ i }, where

W i = w( θ i ) K j=1 w( θ j ) , i = 1, . . . , K,
whose empirical distribution targets the posterior distribution in the sense that

K i=1 W i ϕ( θ i ) → E π(θ|x) (ϕ), for N → ∞,
for any measurable application ϕ. The main challenge when using importance sampling is to find a suitable importance distribution g. If the purpose of g is to be substitute for π(•|x) then the Effective Sample Size (ESS) of the particle cloud must be high enough. The ESS is an indicator taking values between 1 and N that measures the degeneracy of the cloud of particles. It corresponds to the size of an iid sample that would match the empirical variance of the cloud of weighted particles {(W i , θ i ), i = 1, . . . , N }. The ESS is estimated by

ESS ≈ 1 K i=1 W 2 i ,
as suggested in Kong et al. [START_REF] Kong | Sequential imputations and bayesian missing data problems[END_REF].

The sequential Monte Carlo algorithm presented in the next section bypasses the choice of a proposal distribution by constructing a sequence of intermediary distributions while maintaining an appropriate effective sample size.

Sequential Monte Carlo algorithmic details

A sequential Monte Carlo algorithm builds a sequence of distribution π s (θ|x), s = 0, . . . , t starting from the prior distribution π 0 (θ|x) = π(θ) and ending on the posterior π t (θ|x) = π(θ|x). To build this sequence of distribution sequence π s (θ|x), s = 0, . . . , t, we consider simulated annealing, see Neal [START_REF] Neal | Annealed importance sampling[END_REF]. This technique gradually activates the likelihood function as

π s (θ|x) = L(x|θ) τ s π(θ) Z s , s = 0, . . . , t, (16)
where τ s , s = 0, . . . , t is a sequence of real numbers such that 0 = τ 0 < τ 1 < . . . < τ t = 1, and the normalizing constant is given by

Z s = Θ L(x|θ) τ s π(θ)dθ.
The smc algorithm initializes a cloud of particles using the prior distribution as

θ (0) i i.i.d.
∼ π(θ), and

W (0) i = 1 K , for i = 1, . . . , K.
To move from one intermediary distribution π s to the next π s+1 , the smc algorithm takes the cloud of particles {(W s i , θ s i ), i = 1, . . . , K} and apply three operations to get {(W s+1 i , θ s+1 i ), i = 1, . . . , K}.

1. (Reweighting step) This step prepares the current cloud to target the next distribution. A particle θ i s is reweighted by

W s+1 i ∝ w s+1 i = π s+1 (θ s i ) π s (θ s i )
, for i = 1, . . . , K,

where ∝ stands for "proportional to" and the w s+1 i 's are unnormalized weights, useful to estimate the normalizing constant as we shall see later. Because the weights W s+1 1 , . . . , W s+1 K are actually importance weights, the targeted distribution π s+1 is chosen so that the weights satisfy

ESS ≈ 1 K i=1 (W s+1 i ) 2 ≥ ρN ,
where ρ ∈ (0, 1). The selection of the next target reduces to picking a suitable temperature τ s+1 . This is done via binary search and ρ is set to 1/2 following up on the recommendation of Jasra et al. [START_REF] Jasra | Inference for Lévy-driven stochastic volatility models via adaptive sequential monte carlo[END_REF].

2. (Resampling step) Particles θ s 1 , . . . , θ s K are sampled from the particle clouds {(W s+1 i , θ s i ), i = 1, . . . , K}. A simple multinomial resampling is used here, but note that alternative schemes discussed for instance in the work of Gerber et al. [START_REF] Gerber | Negative association, ordering and convergence of resampling methods[END_REF] are also possible.

3. (Move step) Metropolis-Hasting within Gibbs moves are applied to the particles θ s 1 , . . . , θ s K to yield the new generation of particles θ s+1 1 , . . . , θ s+1 K . The matrix H of the Markov Kernel K is given by Σ

• 2.38/ √ d,
where Σ is the empirical variance-covariance matrix of the particles system {(W s+1 i , θ s i ), i = 1, . . . , K}. The number of transitions k ∈ N to be applied is set to ensure the diversification of the particle cloud. In practice, the Markov kernel is applied once to each particle. The acceptance rate p a is estimated after this pilot run and k is then given by

k = max k max , min k min , log(1 -c) log(1 -p a ) ,
where k min and k max denotes the minimum and maximum number of transitions, and c ∈ (0, 1) is the probability that each particle is moved at least once. Note that k min , k max and c are the user-defined parameters of the smc algorithm. The new particles θ s+1 1 , . . . , θ s+1 K are sampled from π s+1 and are equally weighted with

W s+1 i = 1/K for i = 1, . . . , K.
The adaptative choice of the target distribution in step 1 and the calibration of H and k in step 3 are standard smc algorithmic tricks used for instance in the paper of South et al. [START_REF] South | Sequential monte carlo samplers with independent markov chain monte carlo proposals[END_REF] and the smc sampler of the Python package pymc of Salvatier et al. [START_REF] Salvatier | Probabilistic programming in python using PyMC3[END_REF]. The Metropolis-Hasting move step, also used in Nguyen et al. [START_REF] Le | Efficient sequential monte-carlo samplers for bayesian inference[END_REF], is easy to paralellize to optimize the computing time.

The unnormalized weights {w s i , 1 ≤ i ≤ K, 1 ≤ s ≤ t} yield an approximation of the normalizing constant as

Z(x) = Z t = t s=1 Z s Z s-1 ≈ t s=1        1 K K i=1 w i s        .
The accuracy of the estimator depends on the population size (the higher the better) which is chosen by the user according to a computing time budget. The bayes-splicing package implements the sequential Monte Carlo sampler to fit and compare the composite models introduced in Section 2.1. In addition to the posterior probability, it is possible to compute two Bayesian information criteria. The Deviance Information Criterion (DIC), introduced in the work of Spiegelhalter et al. [START_REF] Spiegelhalter | Bayesian measures of model complexity and fit[END_REF] and the Widely applicable Information Criterion (WAIC) of Watanabe [START_REF] Watanabe | Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory[END_REF]. For a comprehensive discussion about the information criteria used in Bayesian statistics, I refer the reader to the work of Gelman et al. [START_REF] Gelman | Understanding predictive information criteria for bayesian models[END_REF]. Although this work focuses on the posterior model probability to compare the competing theories, note that the bayes-splicing includes python methods to compute the DIC and WAIC.

Simulation study

The experiment involves generating claim data from a known composite model and fitting it using multiple composite models. The pdf of the composite models in competition are continuous at the threshold, and combine the Weibull or the Inverse-Weibull distribution for the body to the Inverse-Weibull, Lomax or Log-Logistic distribution for the tail. This makes six models to choose from. Note that the choice of the datagenerating model and competing models is arbitrary. We limit ourselves to 6 models so as not to unreasonably increase the calculation times and for the sake of readability of the graphs on which the results are presented.

The sample size of the artificial data ranges from 500 to 5, 000. For each of the 1, 000 simulation runs, the posterior probability of each model is computed. A first objective is to assess the finite sample consistency of the model selection procedure. Once the composite models have been calibrated, an estimation of the quantiles of the first-line insurer surplus distribution, defined in Section 2.2, is produced and compared to the true value. For the latter application, The average number of claims is λ = 250. The xol priority and limit are set to be the 90% and 99% quantiles of the splicing distribution that generated the data. The quantile estimates resulting from the best model (according to the posterior probabilities) and from the Bayesian model averaging procedure are also compared to the true value. Two cases are considered. In Section 4.1, the model is well-specified, meaning that the data generating model belongs to the set of concurrent models. In Section 4.2, the model is misspecified, which means that the data generating model does not belongs to the set of concurrent models.

When the model is well-specified

The losses are drawn from a Weibull(k = 1/2, β = 1) -Lomax(α = 2.5, σ = 1.5) continuous composite model with a threshold γ = 1.5. The prior assumptions of the competing composite models are given in Appendix B, see Table 3. The prior assumption over the threshold parameter is given by γ ∼ Gamma(1, 1). The number of particles is N = 8, 000. The posterior model probabilities of the competing composite models computed for each of the 1, 000 simulation runs for samples of sizes 500, 1000, 2000, and 5000 are shown in Figure 1.

The posterior probabilities discriminate the composite models having the Inverse-Weibull distribution as

When the model is misspecified

The losses are drawn from a Exp(δ = 1/2) -Burr(α = 1.8, β = 2, σ = 3) continuous composite model with a threshold γ = 2.5. The prior assumptions of the competing composite models are given in Table 3. The number of particles is N = 8, 000. The posterior model probabilities of the competing composite models computed for each of the 1, 000 simulation runs for samples of sizes 500, 1000, 2000, and 5000 are shown in Figure 3.

The highest posterior probabilities are reported for the composite models having the Weibull distribution as body, and more specifically for the composite models having the Log-logistic and Lomax distributions as tail.

This tendency is reinforced when the sample size increases. This result was expected, since the exponential distribution is a special case of the Weibull distribution on the one hand and the Lomax and Log-logistic distributions are special cases of the Burr distribution on the other hand. Figure 4 displays the boxplots of the 1% and 0.5% quantiles of the first line insurer surplus Π D -D depending on the composite model used.

When the model is misspecified then the "Best" and "BMA" models provide the most accurate estimation of the quantiles. Note that the BMA approach also seems to decrease the variance of the estimation. This validates the model selection procedure proposed in this work since model misspecification tends to be the rule in real

Application to the Danish fire insurance data

The Danish fire insurance data are highly valued by actuarial researchers for comparing statistical methods dealing with extreme values. The data was retrieved from the R package SMPracticals of the book by Davison [START_REF] Davison | Statistical Models[END_REF]. The aim of this study is to find the most suitable composite model for these data among all the possible combinations of distribution for the body and the tail, the different types of composite model including "continuous", "discontinuous" and "simple", and the threshold selection method. The composite models are compared based on log marginal likelihood and a quantile-based distance, which we now present. Denote by F(x, θ) the cdf of the composite model and let Q(y, θ) = inf{x ∈ R ; F(x, θ) ≥ y} be its generalized inverse function that we refer to as the quantile function. The distance

W (F, x) = 1 n n i=1 Q i/n, θ -x i:n , (17)
aims at comparing the quantiles estimated via the composite model Q i/n, θ , where θ denotes the mean a posteriori, to the empirical quantiles x i:n for i = 1, . . . , n. The notation W comes from the fact that [START_REF] Gelman | Understanding predictive information criteria for bayesian models[END_REF] corresponds to the empirical counterpart of a variation of the Wasserstein distance between the empirical measure and the probability measure with density f (x, θ) with respect to the Lebesgue measure. The distance [START_REF] Gelman | Understanding predictive information criteria for bayesian models[END_REF] summarizes in a single number the information contained in a Q-Q plot. Information criteria based on the likelihood function, such as log marginal likelihood, often tend to focus on the part of the distribution that concentrates many data points. A distance based on quantiles, like [START_REF] Gelman | Understanding predictive information criteria for bayesian models[END_REF], is more sensitive to extreme values and therefore puts more emphasis on a potential mismatch in the tail of the distribution. In Section 5.1, the threshold parameter is treated like any of the other model parameters. In addition to choosing the most suitable composite model, we want to highlight the impact of the selected composite model on the value of the extreme value threshold within this free-threshold approach. Section 5.2 completes the results of Section 5.1 by considering the case where the value of the extreme value threshold is fixed before estimating the other parameters of the composite model.

Simultaneous estimation of the extreme value threshold

In this section, the threshold parameter of the composite model is estimated together with the other parameters.

The prior assumptions over the threshold and mixing parameters (for the discontinuous composite models) read as follows

γ ∼ Uniform([min(x), max(x)]), and p ∼ Uniform([0, 1]).
The number of particles is N = 10, 000. Figure 5 gives the boxplots of the log marginal likelihood and the quantile distance depending on the type of composite model considered. The discontinuity at the threshold improves the overall fit of the model to the data when examining the marginal likelihood values. The difference in terms of matching the empirical quantiles of the data, however, is less obvious. Figure 6 shows the posterior probabilities, denoted by w which stands for weight, of the various composite models grouped by type. The "discontinuous" and "simple" composite models agree on a threshold of 0.82 as if a discontinuity were indeed present in the data. For the "continuous" composite models, the results are more contrasted even if the threshold takes on fairly low values, oscillating between 0.9 and 1. Let us take as an example the Weibull-Log logistic composite model. Figure 8 shows the posterior distributions of the parameters in the "continuous" and "discontinuous" setting. The shape of the posterior distribution of the threshold in the "discontinuous" case Overall, many threshold values reported in Figure 4 are too low because 95% to 99% of the claim amounts are considered extreme. In my humble opinion, one possibility is to define a balanced composite model with a few parameters, for example, the exponential distribution as the body and the Pareto distribution as the tail. The resulting threshold is 1.83, which is reasonable since 60% of claims are below this threshold. This recalls the procedure of detecting a break in the Q-Q plot of the exponential (see Figure 9a) and Pareto distribution (see Figure 9b). The resulting fit seems fair, especially in the tail, see Figure 9c.

(PSLULFDO 7KHRUHWLFDO the danish fire loss data using the posterior mean.

Separate estimation of the extreme value threshold

Many methods for selecting the extreme value threshold have been proposed in the statistical literature and they lead to very different threshold values as shown in Table 2. The value reported in Table 2 were obtained using the tea3 R package. The purpose of this work is not to discuss the merits of each of these methods but rather to illustrate the ability of the bayes-splicing package to take a fixed threshold value before inferring the rest of the parameters of the composite model. We therefore use an automatic threshold selection method presented in a recent work of Bladt et al. of composite models. This can be explained by the flexibility of the Burr distribution thanks to its three parameters. This result justifies the use of a nonparametric assumption on the distribution of the body of the composite model when the threshold is chosen before the estimation of the other parameters. This makes it possible in particular to overcome the uncertainty surrounding the choice of the distribution of the body of the composite model. We further compare the goodness-of-fit of the composite models with fixed threshold to that of the composite models for which the threshold was set by the algorithm. Figure 11 shows the boxplot of the log marginal likelihood and the quantile distance for the composite models having a Pareto or a generalized Pareto tail. The examination of the marginal likelihood indicates that it is preferable to leave the threshold free while fixing the threshold beforehand improves the distance to the empirical quantiles. This last remark does not seem to hold when the composite model is continuous at the threshold. We conclude that for a "continuous" composite model, performing simultaneous estimation is the way to go because the overall fit is better and the quantile matching is as good as with a prefixed threshold. If the threshold is fixed beforehand, it is preferable to opt for a "discontinuous" composite model because the loss of the degree of freedom associated with the mixing parameter p seriously affects the fit to the data.

Conclusions and perspectives

This paper presents an implementation of a smc sampler to fit and compare composite models within a Bayesian framework. The python package can be installed from pip, see bayes-splicing and the notebooks to reproduce the results can be freely downloaded from the following github repository SMCCompoMo. The Bayesian approach, compared to the frequentist approach, takes into account the uncertainty around the parameter estimates and enables to encapsulate expert knowledge in the prior distribution. smc samplers have three advantages over the standard mcmc algorithm: (1) It avoids the fine tuning of some hyperparameters, (2) it provides an appoximation of the normalizing constant as a byproduct, and (3) it is very easy to paralellize to take advantage of the multi-core processors that equip modern computers. The simulation experiment showed the capacity of the algorithm to identify the model that generated the data. The advantage of using the posterior model probabilities to weight the estimates of the quantile of the insurer's surplus returned by the various models has been demonstrated, in particular in the case where the model generating the data does not belong to the set of concurrent models. The analysis of the Danish fire insurance data highlighted the advantages and disadvantages of taking a discontinuous or continuous composite model, and also of letting the algorithm select a threshold or fix it beforehand.

There are many avenues for future research. The Monte Carlo sequential sampler could be improved by considering a more sophisticated particle perturbation scheme than the random walk one, similar to what is proposed by South et al. [START_REF] South | Sequential monte carlo samplers with independent markov chain monte carlo proposals[END_REF]. It is also possible to combine the estimates of the different models using weights based on cross-validation procedures to improve the predictive power of the Bayesian averaging procedure, see the work of Yao et al. [START_REF] Yao | Using stacking to average bayesian predictive distributions (with discussion)[END_REF]. When the threshold is set before the other parameters, then a flexible model for the body of the composite model and a discontinuity at the threshold can be optimal. Along these lines, a nonparametric assumption could be made over the body distribution. Tancredi et al. [START_REF] Tancredi | Accounting for threshold uncertainty in extreme value estimation[END_REF] used a mixture of uniform distributions, MacDonald et al. [START_REF] Macdonald | A flexible extreme value mixture model[END_REF] used a kernel density estimator, Cabras and Castellanos [START_REF] Cabras | A bayesian approach for estimating extreme quantiles under a semiparametric mixture model[END_REF] used an orthogonal polynomial expansion, and Reynkens et al. [START_REF] Reynkens | Modelling censored losses using splicing: A global fit strategy with mixed erlang and extreme value distributions[END_REF] used a mixture of Erlang distributions.

The question is which of these methods best suits the Bayesian framework proposed here. In view of the importance given to the tail of the loss distribution, an inference method that does not rely on the likelihood function could be considered. Minimum distance estimators that minimize a measure of the deviation between the empirical quantiles and those of the model could be a valid direction. The work of Bernton et al. [START_REF] Bernton | On parameter estimation with the wasserstein distance[END_REF] focuses on parameter estimates that minimize the Wasserstein distance which reduces in our case (iid and univariate data) to a distance between quantiles. Posterior distributions may be obtained by applying an Approximate Bayesian Computation (ABC) algorithm. ABC combined to the Wasserstein distance have been considered in the work of Bernton et al. [START_REF] Bernton | Approximate bayesian computation with the wasserstein distance[END_REF] and applied to aggregated insurance data in the work of Goffard and Laub [START_REF] Goffard | Approximate bayesian computations to fit and compare insurance loss models[END_REF].

A Summary of the Sequential Monte Carlo algorithm Algorithm 1 smc sampler for π(θ|x)

1: Set ρ ∈ (0, 1); k min ∈ N; k max ∈ N; c ∈ (0, 1) 2: Initialize s ← 0 ; π 0 (θ) ← π(θ); for i = 1 → K do 10:

Sample θ i ∼ {θ for i = 1 → K do for i = 1 → K do

B A priori assumptions

Table 3 provides an overview of the prior assumtions used to do the Bayesian fit of the composite models in Sections 4 and 5. Note that the bayes-splicing allows the user to choose different hyper-parametrization than what is specified in Table 3. Note also that beta and uniform distribution can be used as in addition to the

Example 1 .

 1 bayes-splicing package implements three composite models with continuity and differentiability at the threshold including the lognormal-Pareto model, the Weibull-Pareto model and the gamma-Pareto model described in Example 1. If f 1 is the pdf of the gamma distribution Gamma(r, m) and f 2 is the pdf of the Pareto distribution Pareto(α, γ), then continuity and diffrentiability of (1) lead to express m and p in terms of the other parameters as

Figure 1 :WFigure 2 :

 12 Figure 1: Posterior probabilities assigned to each composite model for sample of sizes 500, 1000, 2000, and 5000 drawn from a continuous Weibull -Lomax composite model.

Figure 3 :Figure 4 :

 34 Figure 3: Posterior probabilities assigned to each composite model for sample of sizes 500, 1000, 2000, and 5000 drawn from a continuous Exp -Burr composite model.

Figure 5 :Figure 6 :

 56 Figure 5: Log marginal likelihood and quantile distance of the fitted composite models depending on their type.

Figure 7 :

 7 Figure 7: Posterior mean of the threshold resulting from the simultaneous estimation of the composite model parameters.

Figure 8 :

 8 Figure 8: Posterior plots of the composite model parameters.

Figure 9 :

 9 Figure 9: Q-Q plots of the Exp model, Pareto model, and the Exp-Pareto continuous composite model fitted to

  [START_REF] Bladt | Threshold selection and trimming in extremes[END_REF]. It consists in minimizing the asymptotic mean squared error of a trimmed version of the Hill estimator of the tail index parameter of the Pareto distribution. The procedure leads to a threshold of 2.78. Since the threshold selection method assumes that the tail is of Pareto type, only the composite models having a Pareto or generalized Pareto tail are considered. The number of particles is N = 10, 000. Figure10shows the posterior probabilities, denoted by w which stands for weight, of the various composite models grouped by type.The highest posterior probability is obtained by the composite Burr-Pareto model and this for all types

Figure 10 :

 10 Figure 10: Posterior probabilities of the composite models depending on their type with the fixed threshold approach.

Figure 11 :

 11 Figure11: Log marginal likelihood and quantile distance of the composite models according to their type, and if the threshold was estimated before or at the same time as the other parameters.

3 :=

 3 for i = 1 → K do 4: θ 0 i ∼ π(θ) ; W 0 i ← 1/K 5: end for 6: while π s (θ) π(θ|x) do π s+1 (θ s i )/π s (θ s i ), i = 1, . . . , K 8: Compute Σ = Cov {(W s+1 i , θ s i ), i = 1, . . . , K} 9:

K

  } with probabilities W s+1 j , for 1 ≤ j ≤ K

15 :

 15 ← K H ( θ i , •) where K H ( θ i , •) where H = 2.38 √ Compute p a = N -1 K i=1 I θ * i = θ i ; k = max k max , min k min , log(1-c) log(1-p a )16:

  i , •) where K * (k-1) H ( θ * i , •) corresponds to k -1 Metropolis-Hasting-Gibbs moves 18:

Table 1 :

 1 List of distribution for bulk and the tail of the composite models.

	many previous works in the statistical and actuarial science litterature, the tail component of the composite
	model may not belong to the Pareto distribution family. Three composite model settings are considered. The
	"discontinuous" composite model assumes that the mixing parameter p has a prior distribution (uniform or

Table 2 :

 2 Threshold values depending on the extreme value technique used.

		7\SH GLVFRQWLQXRXV	7\SH VLPSOH	7\SH FRQWLQXRXV
	7DLO	3DUHWR7DLO *3'7DLO		Z
		/RJ/RJLVWLF /RPD[ %XUU ([S *DPPD :HLEXOO ,QYHUVH*DXVVLDQ /RJQRUPDO ,QYHUVH:HLEXOO %RG\	/RJ/RJLVWLF /RPD[ %XUU ([S *DPPD :HLEXOO ,QYHUVH*DXVVLDQ /RJQRUPDO ,QYHUVH:HLEXOO %RG\	/RJ/RJLVWLF /RPD[ %XUU ([S *DPPD :HLEXOO ,QYHUVH*DXVVLDQ /RJQRUPDO ,QYHUVH:HLEXOO %RG\

https://pypi.org/project/bayes-splicing/

scipy.integrate method from the scipy python library https://docs.scipy.org/doc/scipy/reference/generated/scipy. integrate.quad.html

tea package documentation

Acknowledgmements

The author thanks the associate editor and the two anonymous reviewers for their insightful comments which greatly improved the manuscript. The author's work is conducted within the Research Chair DIALog under the aegis of the Risk Foundation, an initiative by CNP Assurances.

Inverse-Gaussian Inverse-Gaussian(µ, λ) µ, λ > 0 µ, λ ∼ Gamma(1, 1)

Inverse-Gamma Inverse-Gamma(r, m) r, m > 0 r, m ∼ Gamma(1, 1)

Table 3: Prior assumptions over the parameters of distribution for the bulk and tail of the composite models.

gamma and normal distributions. The γ parameter in the Pareto and Generalized Pareto cases corresponds to a threshold. This parameter should be set to the minimum of the data points if the models are considered alone.

The prior specification of the mixing parameter p and the threshold parameter γ of the composite models are given in the main text.

C Computing time analysis

This section gives an overview of the time required to fit composite models using the home-made SMC algorithm of the bayes-splicing package. Calculation times are provides for composite models of varying complexity and with changing regularity assumptions at the threshold. The impact of the size of the particle population is studied as well as whether the move step is parallelized or not. The models were fitted to the danish fire insurance dataset. The experiment has been conducted on a Dell Latitude 7290 equipped with an Intel(R) Core(TM) i7-8650U CPU @ 1.90GHz. When the paralellized version of the algorithm is used, 4 cores are mobilized simultaneously. Table 4 shows the computation times associated with fitting discontinuous composite models using the fixed-threshold approach. It is observed that the computation time increases with the number of parameters and the size of the particle population. The parallelization of the move step improves the calculation times, its interest is all the more clear as the number of particles increases. Table 5 shows the computation times associated with fitting continuous composite models using the fixed-threshold approach.

Paralell Model K = 5, 000 K = 10, 000 K = 20, 000 K = 5, 000 K = 10, 000 K = 20, 000 depending on the composite model, the number of particles and whether the move step is paralellized.

Compared to Table 4, the calculation times are slightly lower since one less parameter is to be estimated. Table 6 and Table 7 shows the computation times associated with fitting discontinuous and continuous composite models using the free-threshold approach. Compared to Table 4 and Table 5, the calculation times are slightly longer since one more parameter has to be estimated. Otherwise, the same remarks apply for the impact of the number of parameters and the size of the particle population.

Paralell Model K = 5, 000 K = 10, 000 K = 20, 000 K = 5, 000 K = 10, 000 K = 20, 000 depending on the composite model, the number of particles and whether the move step is paralellized.