GEMAct: a Python package for non-life (re)insurance modeling

This paper introduces $\textbf{gemact}$, a $\textbf{Python}$ package for actuarial modelling based on the collective risk model. The library supports applications to risk costing and risk transfer, loss aggregation, and loss reserving. We add new probability distributions to those available in $\textbf{scipy}$, including the (a, b, 0) and (a, b, 1) discrete distributions, copulas of the Archimedean family, the Gaussian, the Student t and the Fundamental copulas. We provide an implementation of the AEP algorithm for calculating the cumulative distribution function of the sum of dependent, non-negative random variables, given their dependency structure specified with a copula. The theoretical framework is introduced at the beginning of each section to give the reader with a sufficient understanding of the underlying actuarial models.


Introduction
In non-life insurance, the accurate representation and quantification of future losses is a foundational task, central to several areas ranging from pricing and reserving to risk management.Indeed, the actuarial literature is rich in models that are relevant in such applications.Among those, the collective risk model has been widely studied as it is mathematically tractable, it requires little and general information, and it can be efficiently implemented (Klugman, Panjer, & Willmot, 2012;Embrechts & Frei, 2009;Parodi, 2014).In particular, by knowing the frequency and severity distributions of the losses, it is possible to compute the distribution of the aggregate (total) losses.The collective risk model represents the common thread of this work and we developed gemact to provide a collection of tools for (re)insurance modelling under a unified formal framework.
After a brief discussion on how to install the software in Section 2, we introduce the statistical framework of the collective risk model in Section 3. There, we define an aggregate loss distribution as a random sum of i.i.d.random variables, which can be computed using the recursive formula (Panjer, 1981), the discrete Fourier transform (Bühlmann, 1984;S. Wang, 1998;Grübel & Hermesmeier, 1999) and a Monte Carlo simulation approach (Klugman et al., 2012, p. 467).Once the aggregate loss distribution is available, its expected value can be used for costing purposes.In this respect, the package supports typical coverage modifiers like (individual and aggregate) deductibles, limits, and reinstatements (see Sundt, 1990).Also, we consider different methods for the discretisation of continuous distributions (Gerber, 1982).
Often, it is necessary to model the sum of a fixed number of dependent random variables.In order to do so, in Section 4, we introduce the AEP algorithm (Arbenz, Embrechts, & Puccetti, 2011) and a Monte Carlo simulation approach for evaluating the cumulative distribution function of a sum of random variables with a given dependency structure.The dependency structure can be specified with the copulas we implemented.These are listed in Appendix A and include copulas of the Archimedean family, the Gaussian, the Student t and the Fundamental copulas (Nelsen, 2007).
Lastly, assuming a collective risk model holds for the cells of a loss development triangle, it is possible to define the stochastic claims reserving model in Ricotta and Clemente (2016); Clemente, Savelli, and Zappa (2019).In this case, the user obtains information on the frequency and severity parameters of the cells from the Fisher-Lange method (Fisher & Lange, 1973).Both these approaches are described in Section 5.

Context, scope, and contributions
In the recent years, programming languages and statistical computing environments, such as Python (Van Rossum & Drake, 2009) and R (R Core Team, 2017), have become increasingly popular (Ozgur, Colliau, Rogers, & Hughes, 2022).Currently, coding skills form part of the body of knowledge of actuaries and actuarial science researchers.In R, an extensive implementation for aggregate loss modelling based on the collective risk theory is the actuar package (Dutang, Goulet, & Pigeon, 2008;Dutang, Goulet, & Langevin, 2022).An available library in Python is aggregate, which implements the computation of compound probability distributions via fast Fourier transform convolution algorithm (Mildenhall, 2022).This package employs a specific grammar for the user to define insurance policy declarations and distribution calculation features.Direct access to its objects and their components is also possible.
With regards to claims reserving, chainladder offers in Python standard aggregate reserving techniques, like deterministic and stochastic chain-ladder methods, the Bornhuetter-Ferguson model, and the Cape Cod model (Bogaardt, 2022).This package is available in R and Python (Gesmann et al., 2022).Furthermore, apc provides the family of age-period-cohort approaches for reserving.This package is also available in both the above-mentioned programming languages (Nielsen, 2015).
When it comes to the topic of dependence modelling via copulas, in Python one can use the copulas and copulae packages (Bok, 2022;Lab, 2022).Similarly, copula features in R are implemented in copula, see the package and its extensions in Jun Yan (2007); Ivan Kojadinovic and Jun Yan (2010); Marius Hofert and Martin Mächler (2011).
In this manuscript, we present an open-source Python package that extends the existing software tools available to the actuarial community.Our work is primarily aimed at the academic audience, who can benefit from our implementation for research and teaching purposes.Nonetheless, our package can also support non-life actuarial professionals in prototypes modelling, benchmarking and comparative analyses, and ad-hoc business studies.
From the perspective of the package design, gemact adopts an explicit, direct and integrated object-oriented programming (OOP) paradigm.In summary, our goal is to provide: • A computational software environment that gives users control over mathematical aspects and actuarial features, enabling the creation of models tailored to specific needs and requirements.
• An object oriented system whose elements (i.e.objects, methods, attributes) can be accessed and managed via our API in such a way as to be interactive and suitable for users familiar with OOP designs and with the underlying modelling framework.
• A collection of extensible libraries to make gemact survive over time.Our package is designed in an attempt to be easily extended or integrated with new functionalities, modules, and to respect the attributes that qualify extensible software (Johansson & Löfgren, 2009).Namely, a modification to the functionalities should involve the least number of changes to the least number of possible elements (modifiability, Bass, Clements, & Kazman, 2003, p. 137), the addition of new requirements should not raise new errors (maintainability, Sommerville, 2011, p. 25), and the system should be able to expand in a chosen dimension without major modifications to its architecture (scalability, Bondi, 2000).
From the perspective of actuarial advancements and developments, gemact provides an implementation to established algorithms and methodologies.In particular, our package: • Implements (a, b, 0) and (a, b, 1) classes of distributions for describing the loss frequency (Klugman et al., 2012, p. 505), and further continuous distributions to model the loss severity, like the generalised beta (Klugman et al., 2012, p. 493).Additional details can be found in Appendix A. Moreover, it integrates these into scipy distributions (Virtanen et al., 2020).• Offers the first open-source software implementation of the AEP algorithm (Arbenz et al., 2011) for evaluating the cumulative distribution function of a sum of random variables with a given dependency structure specified via a copula.• Includes the Student t copula and a method for numerically approximating its cumulative distribution function (Genz & Bretz, 1999, 2002).• Implements the stochastic claims reserving model described by Ricotta and Clemente (2016); Clemente et al. (2019) based upon the collective risk model apparatus.

Installation
The production version of the package is available on the Python Package Index (PyPi).Users can install gemact via pip, using the following command in the operating system command line interface.
Examples on how to get started and utilise objects and functionalities of our package will be shown below.This work refers to production version 1.2.1.

Loss model
Within the framework of the collective risk model (Embrechts & Frei, 2009), all random variables are defined on some fixed probability space (Ω, F, P ).Let • N be a random variable taking values in N 0 representing the claim frequency.
• {Z i } i∈N be a sequence of i.i.d non-negative random variables independent of N ; Z is the random variable representing the individual (claim) loss.
The aggregate loss X, also referred to as aggregate claim cost, is with 0 i=1 Z i = 0. Details on the distribution functions of N , Z and X are discussed in the next sections.Equation ( 1) is often referred to as the frequency-severity loss model representation.This can encompass common coverage modifiers present in (re)insurance contracts (Parodi, 2014, p. 50).More specifically, let us consider: • For a ∈ [0, 1], the function Q a apportioning the aggregate loss amount: (2) • For c, d ≥ 0, the function L c,d applied to the individual claim loss: Herein, for each and every loss, the excess to a deductible d (sometimes referred to as priority) is considered up to a cover or limit c.In line with Albrecher, Beirlant, and Teugels (2017, p. 34), we denote [d, d + c] as layer.An analogous notation is found in Ladoucette and Teugels (2006) and Parodi (2014, p. 46).Similarly to the individual loss Z i , Equation (3) can be applied to the aggregate loss X.
Computing the aggregate loss distribution is relevant for several actuarial applications (Parodi, 2014, p. 93).
The gemact package provides various methods for calculating the distribution of the loss model in Equation ( 1) that allow the transformations of Equation (2), Equation (3) and their combinations to be included.

Risk costing
In this section, we describe an application of the collective risk model of Equation ( 1).The expected value of the aggregate loss of a portfolio constitutes the building block of an insurance tariff.This expected amount is called pure premium or loss cost and its calculation is referred as risk costing (Parodi, 2014, p. 282).Insurers frequently cede parts of their losses to reinsurers, and risk costing takes this transfers into account.Listed below are some examples of basic reinsurance contracts whose pure premium can be computed using gemact.
• The Quota Share (QS), where a share a of the aggregate loss ceded to the reinsurance (along with the respective premium) and the remaining part is retained: • The Excess-of-loss (XL), where the insurer cedes to the reinsurer each and every loss exceeding a deductible d, up to an agreed limit or cover c, with c, d ≥ 0: (5) • The Stop Loss (SL), where the reinsurer covers the aggregate loss exceedance of a (aggregate) deductible v, up to a (aggregate) limit or cover u, with u, v ≥ 0: The model introduced by Equation (1) and implemented in gemact can be used for costing contracts like the Excess-of-loss with Reinstatements (RS) in Sundt (1990).Assuming the aggregate cover u is equal to (K + 1)c, with K ∈ Z + : where K is the number of reinstatement layers and l k ∈ [0, 1] is the reinstatement premium percentage, with k = 1, . . ., K. When l k = 0, the k-th resinstatement is said to be free.In detail, the logic we implemented implies that whenever a layer is used the cedent pays a reinstatement premium, i.e. l k P RS , and the cover c is thus reinstated.The reinstatement premium will usually be paid in proportion to the amount that needs to be reinstated (Parodi, 2014, p. 52).In practice, the reinstatement premium percentage l k is a contractual element, given a priori as a percentage of the premium paid for the initial layer.In fact, in gemact l k is provided by the user.The mathematics behind the derivation of P RS is beyond the scope of this manuscript; the interested reader can refer to Sundt (1990), Parodi (2014, p. 325) and Antal (2009, p. 57).

Computational methods for the aggregate loss distribution
The cumulative distribution function (cdf) of the aggregate loss in Equation ( 1) is: where Moreover, the characteristic function of the aggregate loss ϕ X : R → C can be expressed in the form: where P N (t) = E t N is the probability generating function of the frequency N and ϕ Z (t) is the characteristic function of the severity Z (Klugman et al., 2012, p. 153).
The distribution in Equation ( 8), except in a few cases, cannot be computed analytically and its direct calculation is numerically expensive (Parodi, 2014, p. 239).For this reason, different approaches have been analysed to approximate the distribution function of the aggregate loss including parametric and numerical quasi-exact methods (for a detailed treatment refer to Shevchenko, 2010).Amongst the latter, gemact implements Monte Carlo simulation (Klugman et al., 2012, p. 467), discrete Fourier transform (Bühlmann, 1984;S. Wang, 1998;Grübel & Hermesmeier, 1999) and the so-called recursive formula (Panjer, 1981).A brief comparison of accuracy, flexibility and speed of these methods can be found in Parodi (2014, p. 260) and Wüthrich (2023, p. 127).This section details these last two computational methods based on discrete mathematics to approximate the aggregate loss distribution.

Discrete Fourier transform
The discrete Fourier transform (DFT) of the severity {f } is, for k = 0, ..., m − 1, the sequence where The original sequence can be reconstructed with the inverse DFT: The sequence of probabilities {g} = {g 0 , g 1 , . . ., g m−1 } can be approximated taking the inverse DFT of The original sequence can be computed efficiently with a fast Fourier transform (FFT) algorithm, when m is a power of 2 (Embrechts & Frei, 2009).

Recursive formula
Assume that the frequency distribution belongs to the (a, b, 0) class, i.e. for k ≥ 1 and a, b ∈ R: Here, p 0 is an additional parameter of the distribution (Klugman et al., 2012, p. 505).The (a, b, 0) class can be generalised to the (a, b, 1) class assuming that the recursion in Equation ( 12) holds for k = 2, 3, . . .
The recursive formula was developed to compute the distribution of the aggregate loss when the frequency distribution belongs to the (a, b, 0) or the (a, b, 1) class.The sequence of probabilities {g} can obtained recursively using the following formula: with 1 ≤ s ≤ m − 1 and given the initial condition g 0 = P N (f 0 ).

Severity discretisation
The calculation of the aggregate loss with DFT or with the recursive formula requires an arithmetic severity distribution (Embrechts & Frei, 2009).Conversely, the severity distribution in Equation ( 1) is often calibrated on a continuous support.In general, one needs to choose a discretisation approach to arithmetise the original distribution.This section illustrates the methods for the discretisation of a continuous severity distribution available in the gemact package.
Let F z : R + → [0, 1] be the cdf of the distribution to be discretised.For a bandwidth, or discretisation step, h > 0 and an integer m, a probability f j is assigned to each point hj, with j = 0, . . ., m − 1.Four alternative methods for determining the values for f j are implemented in gemact.
1.The method of mass dispersal: 2. The method of the upper discretisation: 3. The method of the lower discretisation: 4. The method of local moment matching:   The default discretisation method in gemact is the method of the mass dispersal, in which each hj point is assigned with the probability mass of the h-span interval containing it (Figure 1 top left).The upper discretisation and lower discretisation methods generate, respectively, pointwise upper and lower bounds to the true cdf.Hence, these can be used to provide a range where the original F z is contained (see Figure 1, top right and bottom left graphs).The method of local moment matching (Figure 1 bottom right) allows the moments of the original distribution F z to be preserved in the arithmetic distribution.A more general definition of this approach can be found in Gerber (1982).We limited this approach to the first moment as, for higher moments, it is not well defined and it can potentially lead to a negative probability mass on certain lattice points (Embrechts & Frei, 2009).The method of local moment matching shall be preferred for large bandwidths.However, considering that there is no established analytical procedure to determine the optimal step, we remark that the choice of the discretisation step is an educated guess and it really depends on the problem at hand.In general, h should be chosen such that it is neither too small nor too large relative to the severity losses.In the first case, the hj points are not sufficient to capture the full range of the loss amount and the probability in the tail of the distribution exceeding the last discretisation node h(m − 1) is too large.In the second case, the granularity of the severity distribution is not sufficient, and small losses are over-approximated.Additional rules of thumb and guidelines for the choice of discretisation parameters can be found in Parodi (2014, p. 248).For example, one option is to perform the calculation with decreasing values of h and check, graphically or according to a predefined criterion, whether the aggregate distribution changes substantially (Embrechts & Frei, 2009).The reader should refer to Klugman et al. (2012, p. 179) and Embrechts and Frei (2009) for a more detailed treatment and discussion of discretisation methods.
The above-mentioned discretisation methods are modified accordingly to reflect the cases where the transformation of Equation ( 3) is applied to the severity (Klugman et al., 2012, p. 517).Below, we illustrate how to perform severity discretisation in gemact.
Once a continuous distribution is selected from those supported in our package (see Appendix A), the severity distribution is defined via the Severity class.
>>> from gemact.lossmodelimport Severity >>> severity = Severity(dist='gamma', par={'a': 5}) The dist argument contains the name of the distribution, and the par argument specifies, as a dictionary, the distribution parameters.In the latter, each item key-value pair represents a distribution parameter name and its value.Refer to distributions module for a list of the distribution names and their parameter specifications.
The discretize method of the Severity class produces the discrete severity probability sequence according to the approaches described above.Below, we provide an example for mass dispersal.
• The number of nodes (m) set in the n_discr_nodes argument.
• The severity discretisation step (h) is in the discr_step argument.
• If necessary, a deductible specifying where the discretisation begins.The default value is zero.
After the discretisation is achieved, the mean of the discretised distribution can be calculated.
>>> import numpy as np >>> discrete_mean = np.sum(massdispersal['nodes']* massdispersal['fj']) >>> print('Mean of the discretised distribution:', discrete_mean) Mean of the discretised distribution: 5.000000000000079 Additionally, the arithmetic distribution obtained via the severity discretisation can be visually examined using the plot_discretized_severity_cdf.This method is based on the pyplot interface to matplotlib (Hunter, 2007).Hence, plot_discretized_severity_cdf can be used together with pyplot functions, and can receive pyplot.plot arguments to change its output.In the following code blocks we adopt the plot_discretized_severity_cdf method in conjunction with the plot function from matplotlib.pyplot to compare the cdf of a gamma distribution with mean and variance equal to 5, with the arithmetic distribution obtained with the method of mass dispersal above.We first import the gamma distribution from the distributions module and compute the true cdf.
The output of the previous code is shown in the top left graph of Figure 1.To obtain the other graphs, simply set discr_method to the desired approach, i.e. 'upper_discretisation', 'lower_discretisation' or 'localmoments'.

Supported distributions
The gemact package makes for the first time the (a, b, 0) and (a, b, 1) distribution classes (Klugman et al., 2012, p. 81) available in Python.In the following code block we show how to use our implementation of the zero-truncated Poisson from the distributions module.

>>> ztpois = distributions.ZTPoisson(mu=2)
Each distribution supported in gemact has various methods and can be used in a similar fashion to any scipy distribution.Next, we show how to compute the approximated mean via Monte Carlo simulation, with the random generator method for the ZTPoisson class.
>>> from gemact import copulas >>> gumbel_copula = copulas.GumbelCopula(par=1.2, dim=2) >>> values = np.array([[.5, .5]])>>> print('Gumbel copula cdf: ', gumbel_copula.cdf(values)[0])Gumbel copula cdf: 0.2908208406483879 In the above example, it is noted that the copula parameter and dimension are defined by means of the par and the dim arguments, respectively.The argument of the cdf method must be a numpy array whose dimensions meet the following requirements.Its first component is the number of points where the function shall be evaluated, its second component equals the copula dimension (values.shape of the example is in fact (1,2)).
The complete list of the distributions and copulas supported by gemact is available in Appendix A. We remark that the implementation of some distributions is available in both gemact and scipy.stats.However, the objects of the distributions module include additional methods that are specific to their use in actuarial science.Examples are lev and censored_moment methods, which allow the calculation of the limited expected value and censored moments of continuous distributions.Furthermore, the choice of providing a Severity class is in order to have a dedicated object that includes functionalities relevant only for the calculation of a loss model, and not for distribution modelling in general.An example of this is the discretize method.A similar reasoning applies to the Frequency class.In this case, a LossModel object entails Frequency, Severity and PolicyStructure class instances.These correspond to the frequency model, the severity model and the policy structure, respectively.The latter, in particular, is in turn specified via one or more Layer objects, which include coverage modifiers of each separate policy component.

Illustration lossmodel
The following are examples of how to get started and use lossmodel module and its classes for costing purposes.As an overview, Figure 2 schematises the class diagram of the lossmodel module, highlighting its structure and the dependencies of the LossModel class.
The Frequency and the Severity classes represent, respectively, the frequency and the severity components of a loss model.In these, dist identifies the name of the distribution, and par specifies its parameters as a dictionary, in which each item key-value pair corresponds to a distribution parameter name and value.Please refer to distributions module for the full list of the distribution names and their parameter specifications.
The code block below shows how to initiate a frequency model.
>>> from gemact.lossmodelimport Frequency >>> frequency = Frequency( dist='poisson', par={'mu': 4}, threshold=0 ) In practice, losses are reported only above a certain threshold (the reporting threshold) and the frequency model can be estimated only above another, higher threshold called the analysis threshold (Parodi, 2014, p. 323).This can be specified in Frequency with the optional parameter threshold, whose default value is 0 (i.e. the analysis threshold equals the reporting threshold).Severity models in gemact always refer to the reporting threshold.
A loss model is defined and computed through the LossModel class.Specifically, Frequency and Severity objects are assigned to frequency and severity arguments of LossModel to set the parametric assumptions of the frequency and the severity components.Below, we use the severity object we instanced in Section 3.3.In the previous example, in more detail, lm_mc object adopts the Monte Carlo simulation for the calculation of the aggregate loss distribution.This approach is set via the aggr_loss_dist_method equal to 'mc'.The additional parameters required for the simulation are: • the number of simulations n_sim, • the (pseudo)random number generator initialiser random_state.
The cdf of the aggregate loss distribution can be displayed with the plot_dist_cdf method.Moreover, a recap of the computation specifications can be printed with the print_aggr_loss_specs method.When the use_dist argument is set to True (default value) the quantity is derived from the approximated aggregate loss distribution (dist property).Conversely, when it is False, the calculation relies on the closedform formulas of the moments of the aggregate loss random variable.These can be obtained directly from the closed-form moments of the frequency and the severity model transformed according to the coverage modifiers (Parodi, 2014, p. 322).This option is available for mean, std, var (i.e. the variance), coeff_variation (i.e. the coefficient of variation) and skewness methods.It should be noted that the calculation with use_dist=False is not viable when aggregate coverage modifiers are present.In such circumstance, the method must necessarily be based on the approximated aggregate loss distribution.In any situations, it is possible to get the moments of the approximated aggregate loss distribution via the moment method.The central and n arguments specify, respectively, whether the moment is central and the order of the moment.

>>>
>>> lm_mc.moment(central=False,n=1) 19.963155575580227 Furthermore, for the aggregate loss distribution, the user can simulate random variates via the rvs method.The quantile and the cdf functions can be computed via the ppf and the cdf methods.Below is an example of the ppf method returning the 0.80-and 0.70-level quantiles.
>>> lm_mc.ppf(q=[0.80,0.70]) array ([28.81343738, 24.85497983]) The following code block shows the costing of a '20 excess 5' XL reinsurance contract.It can be noted that, in the previous code block, we determined the aggregate loss distribution with 'fft' as aggr_loss_dist_method.In such case, LossModel requires additional arguments for defining the computation process.Namely: • The number of nodes of the aggregate loss distribution n_aggr_dist_nodes.
• The method to discretise the severity distribution sev_discr_method.Above, we opted for the method of mass dispersal ('massdispersal').• The number of nodes of the discretised severity n_sev_discr_nodes (optional).
• The discretisation step sev_discr_step (optional).When a cover is present, gemact automatically adjusts the discretisation step parameter to have the correct number of nodes in the transformed severity support.
The same arguments shall be specified when computing the aggregate loss distribution with the recursive formula, i.e. aggr_loss_dist_method set to 'recursive'.
The costing specifications of a LossModel object can be accessed with the method print_costing_specs().
>>> The next example illustrates the costing of an XL with reinstatements.The PolicyStructure object is set as follows.
Below, we compute the pure premium of Equation ( 7), given the parametric assumptions on the frequency and the severity of the loss model.A PolicyStructure object can handle multiple layers simultaneously.These can overlap and do not need to be contiguous.The length property indicates the number of layers.The code block below deals with three Layer objects, the first without aggregate coverage modifiers and (participation) share of 0.5, the second with reinstatements, the third with an aggregate cover.

Comparison of the methods for computing the aggregate loss distribution
In this section we analyse accuracy and speed of the computation of the aggregate loss distribution using fast Fourier transform (FFT), recursive formula (recursion), and Monte Carlo simulation (MC) approaches, as the number of nodes, the discretisation step and the number of simulations vary.
For this purpose, a costing example was chosen such that the analytical solutions of the moments of the aggregate loss distributions are known and can be compared to those obtained from the approximated aggregate loss distribution.The values of the parameters of the severity and frequency models are taken from the illustration in Parodi (2014, p. 262).Specifically, these and the policy structure specifications are as follows.
• Severity: lognormal distribution with parameters shape = 1.3 and scale = 36315.49,hence whose mean and standard deviation are 84541.68and 177728.30,respectively.
• Frequency: Poisson distribution with parameter µ = 3, with analysis threshold d.It belongs to the (a, b, 0) family described in Section 3.2.2 with parameters a = 0, b = 3 and p 0 = e −3 .
In particular, the accuracy in the approximation of the aggregate loss distribution has been assessed using the relative error in the estimate of the mean, coefficient of variation (CoV) and skewness, with respect to their reference values (i.e.error = estimate/reference -1).The latter are obtained using the following closed-form expressions (Bean, 2000, p. 382).
Their values for the example are in the top part of Table 1.The test of the speed of our implementation has been carried out by measuring the execution time of the approximation of the aggregate loss distribution function with the built-in timeit library.In line with best practice (Martelli, Ravenscroft, & Ascher, 2005, Chapter 18), the observed minimum execution time of independent repetitions of the function call was adopted.
The results of the analysis are reported in  Table 1: Accuracy and speed of the approximation of the aggregate loss distribution using fast Fourier transform (FFT), the recursive formula (recursion), and the Monte Carlo simulation (MC) when varying the number of nodes (m), the discretisation step (h) and number of simulations.The upper table contains the reference values obtained from the closed-form solutions.The lower table reports the execution times in second and the relative errors with respect to the reference values.(Embrechts & Frei, 2009).Furthermore, for this example, in both FFT and recursion, when h increases, the error reduces.Finally, MC approach lies in between the two other alternatives, when it comes to computing time.

Comparison with aggregate FFT implementation
The aggregate package (Mildenhall, 2022) allows to compute the aggregate loss distribution using FFT.
In this section, we compare our implementation with aggregate implementation to show that the two provide similar results.We adopt the same underlying frequency and severity assumptions of Section 3.6, and contracts with different combinations of individual and aggregate coverage modifiers.In particular, we first consider no reinsurance, then an Excess-of-loss (XL), with individual-only coverage modifiers, a Stop Loss (SL) and finally an Excess-of-loss with individual and aggregate coverage modifiers (XL w/agg.).In line with the example in Parodi (2014, p. 262), individual coverage modifiers are c = 1000000 and d = 10000, aggregate coverage modifiers are u = 1000000 and v = 50000.The number of nodes m is set to 2 22 in all the calculations.
The comparison of the speed of the two implementations has been carried out by means of the built-in timeit library.In particular we measured the execution time of both the initialisation of the main computational object and the calculation of the aggregate loss distribution, in order to make the comparison consistent and adequate.In line with best practice (Martelli et al., 2005, Chapter 18), the observed minimum execution time of independent repetitions of the function call was adopted.
As reported by Table 2, the two implementations generate consistent results; their estimates for mean, CoV and skewness tend to coincide for all contracts and are close to the reference values when these are available.When it comes to computing time, gemact takes a similar time, just under one second, for all the examples considered.For the case without reinsurance and the XL, aggregate performs slightly better.Conversely, for SL and XL w/agg., gemact is more than twice as fast.

Loss aggregation
In insurance and finance, the study of the sum of dependent random variables is a central topic.A notable example is risk management, where the distribution of the sum of certain risks needs to be approximated and analysed for solvency purposes (Wilhelmy, 2010).Another application is the pricing of financial and (re)insurance contracts where the payout depends on the aggregation of two or more dependent outcomes (see for example Cummins, Lewis, & Phillips, 1999; P. Wang, 2013).In this section, in contrast with the collective risk theory in Section 3, we model the sum of a given number of random variables d > 1 that are neither independent nor necessarily identically distributed.
More specifically, consider now the random vector is known analytically or can be numerically evaluated in an efficient way.For a real threshold s, the gemact package implements the AEP algorithm and a Monte Carlo simulation approach to model given a set of parametric assumptions on the one-dimensional marginals X 1 , ..., X d and their copula.
More specifically, the AEP algorithm is designed to approximate Equation ( 14) through a geometric procedure, without relying on simulations or numerical integration of a density.In Appendix B a brief description of the algorithm is given.For a complete mathematical treatment of the subject, the reader should refer to Arbenz et al. (2011).

Illustration lossaggregation
Below are some examples of how to use the LossAggregation class.This belongs to the lossaggregation module, whose class diagram is depicted in Figure 3.The main class is LossAggregation, which is the computation object of the random variable sum.This depends on the classes Margins and Copula.Evidently, the former represents the marginal distributions and the latter describes the dependency structure, i.e. the copula.
Consistent with the gemact framework, the specifications needed to instantiate Margins and Copula are akin to those of the Frequency and Severity classes in lossmodel.Copula objects are specified by: • dist: the copula distribution name as a string (str), • par: the parameters of the copula, as a dictionary.
Likewise, Margins objects are defined by: • dist: the list of marginal distribution names as strings (str), • par: the list of parameters of the marginal distributions, each list item is a dictionary.
Please refer to Table 7 and Table 8 of Appendix A for the complete list of the supported distributions and copulas.
In the following code block, we show how to calculate the cdf of the sum of a generalised Pareto and a lognormal dependent random variables, using the AEP algorithm ('aep') and the Monte Carlo simulation approach ('mc').The underlying dependency structure is a Frank copula.LossAggregation includes other functionalities like the survival function sf, the quantile function ppf and the generator of random variates rvs.Furthermore, for the Monte Carlo simulation approach, it is possible to derive empirical statistics and moments using methods such as for example moment, mean, var, skewness, censored_moment.To conclude, the last code block illustrates the calculation of the quantile function via the ppf method.

Comparison of the methods for computing the cdf
In this section, we compared our implementation of the AEP algorithm with the alternative solution based on Monte Carlo simulation in terms of speed and accuracy.Accuracy in the calculation of the cdf has been assessed by means of the relative error with respect to the reference value for a chosen set of quantiles (i.e.error = estimate/reference -1).We replicate the experiment in Arbenz et al. (2011) and take the reference values that the authors computed in the original manuscript.The analysis considers four different Clayton-Pareto models, for d = 2, 3, 4, 5, with the following parametric assumptions.In the two-dimensional case (d = 2), the tail parameters γ of the marginal distributions are 0.9 and 1.8; the Clayton copula has parameter θ = 1.2.In three dimensions (d = 3), the additional marginal has parameter γ = 2.6, and the copula has θ = 0.4.For the four-dimensional (d = 4) and five-dimensional (d = 5) cases, the extra marginal component has parameter equal to 3.3 and 4, and the Clayton copula has parameter 0.2 and 0.3, respectively.The cdf has been evaluated at the quantiles s = {10 0 , 10 2 , 10 4 , 10 6 } for d = 2 and d = 3, and s = {10 1 , 10 2 , 10 3 , 10 4 }, when d = 4 and d = 5.
The test of the speed of our implementation has been carried out by measuring the execution time of the cdf function for a single quantile with the built-in timeit library.In line with best practice (Martelli et al., 2005, Chapter 18), the observed minimum execution time of independent repetitions of the function call was adopted.
Table 3 shows the results of our comparison of the two methodologies for the calculation of the cdf.
It can be noted that our implementation of the AEP algorithm shows a high accuracy in the calculation of the cdf, for all quantiles and dimensions.Its precision also remains valid for the five-dimensional case.The figures are in line with the results of the original study.In general, in cases considered, the AEP algorithm is closer to the reference values than the Monte Carlo simulation approach.Nevertheless, the latter shows contained errors whose order of magnitude is 10 −2 at most, when the number of simulation is set to the lowest value.On the other hand, the AEP algorithm is outperformed by the Monte Carlo simulation approach, in terms of execution speed.The largest gaps have been observed especially when the number of iterations of the AEP is higher and the dimension is 4 and 5.However, it should be noted that the computational times for the AEP algorithm remains acceptable, in most cases below one second even in high dimensions.For the  3: Accuracy and speed of the cdf calculation using the AEP algorithm and the Monte Carlo simulation approach (MC) for the sum of Pareto random variables coupled with a Clayton copula for different dimensions and quantiles.The upper table contains the reference values from Arbenz et al. (2011).The lower table reports the execution times in seconds and the relative errors with respect to the reference values.The time column on the left is about the execution times of gemact, while the time column on the right (labelled with a *) lists those of the original manuscript.sake of completeness and as a reference, the computational time figures of the first implementation, reported in the study of the original manuscript, are also given.
To conclude, we perform a sensitivity study of the two above-mentioned approaches for different dependency structures and number of dimensions.Figure 4 shows how the cdf value calculated with the AEP algorithm changes as the underlying degree of dependency varies, in the cases of a bivariate Gaussian copula and a three-dimensional Clayton copula.The solid black lines of the graphs represent the cdf evaluated at four quantiles s for increasing values of ρ (the non-diagonal entry of the correlation matrix) and θ parameters, for the Gaussian copula and the Clayton copula respectively.The cdf values were also compared with those obtained by the Monte Carlo simulation approach.The average absolute difference between the results of the two methods, across the four quantiles, is highlighted by the dotted red line.It can be seen that this remains stable at low values, regardless of the underlying dependency structure.In all cases, the results produced by the two methods almost coincide.

|AEP -MC|
Figure 4: Sensitivity analysis of the cdf at four quantiles s calculated using the AEP algorithm (n_iter = 7) for different copula models, dimensionality, and underlying degree of dependency.The values of s are 1.25, 1.5, 1.75, 1.95 for the bivariate Gaussian copula (left plot) and 2, 2.3, 2.65, 2.85 for the three-dimensional Clayton (right plot).Each solid black line indicates the values of the cdf for a given s, as the respective parameters ρ and θ change.The results of the AEP algorithm correspond with those of the Monte Carlo simulation approach (MC), using 10 7 number of simulations.The dashed red line represents the average absolute difference between the two method cdf values, calculated across the four quantiles.

Loss reserve
In non-life insurance, contracts do not settle when insured events occur.At the accident date, the insured event triggers a claim that will generate payments in the future.The task of predicting these liabilities is called claims reserving and it assesses claim outstanding loss liabilities (see for example Wüthrich & Merz, 2015, p. 11).In the present work, we refer to the total outstanding loss liabilities of past claims as the loss reserve or claims reserve. Figure 5 sketches an example of the timeline evolution of an individual claim.The insured event occurs within the insured period but the claim settles after several years.In particular, after the claim is reported, the insurance company makes an initial quantification of the claim payment size, the so-called case estimate.Two payments occur thereafter.These payments are not known at the evaluation date and they require to be estimated.In between the payments, the case estimate is updated.Until a claim settle, the insurance company refers to it as an open claim.In certain circumstances, settled claims can be reopened (Friedland, 2010, p. 431).

Insured period
Cash flows

Accident Date
Reporting Date Settlement In this section, we first define the development triangles, the data structure commonly used by actuarial departments for claims reserving (Friedland, 2010, p. 51).These are aggregate representations of the individual claim data.Then, we present the collective risk model for claims reserving in Ricotta and Clemente (2016) and Clemente et al. (2019), which allows to estimate the variability of the reserve.This model requires extra-inputs from a deterministic model to be implemented.In this manuscript, we rely on the Fisher-Lange model (Fisher & Lange, 1973) as proposed in Savelli and Clemente (2014).Further details can be found in Appendix C.

Problem framework
Let the index i = 0, . . ., J denotes the claim accident period, and let the index j = 0, . . ., J represents the claim development period, over the time horizon J > 0. The so-called development triangle is the set Below is the list of the development triangles used in this section.
• The triangle of incremental paid claims: with x i,j being the total payments from the insurance company for claims occurred at accident period i and paid in period i + j.
• The triangle of the number of paid claims: with n i,j being the number of claim payments occurred in accident period i and paid in period i + j.The triangle of average claim cost can be derived from the incremental paid claims triangle and the number of paid claims.Indeed, we define: where m i,j is the average claim cost for accident period i and development period j.
• The triangle of incremental amounts at reserve: R (T ) = {r i,j : (i, j) ∈ T } , with r i,j being the amount booked in period i + j for claims occurred in accident period i.
• The triangle of the number of open claims: with o i,j being the number of claims that are open in period i + j for claims occurred in accident period i.
• The triangle of the number of reported claims: with d i,j being the claims reported in period j and belonging to accident period i.Often, the number of reported claims is aggregated by accident period i.
We implemented the model of Ricotta and Clemente (2016); Clemente et al. (2019), hereafter referred to as CRMR, which connects claims reserving with aggregate distributions.Indeed, the claims reserve is the sum of the (future) payments The authors represent the incremental payments in each cell of the triangle of incremental payments as a compound mixed Poisson-gamma distribution under the assumptions in (Ricotta & Clemente, 2016).For i, j = 0, . . ., J , where the index h is referred to the individual claim severity.The random variable ψ follows a gamma distribution, see Table 4, and introduces dependence between claim-sizes of different cells.

Predicting the claims reserve
In order to estimate the value for the claims reserve R, this approach requires additional parametric assumptions on N i,j and Z h;i,j when i + j > T .In this section, we illustrate how to use the results from the Fisher-Lange to determine the parameters of N i,j and Z h;i,j .The Fisher-Lange is a deterministic average cost method for claims reserving (Institute and Faculty of Actuaries, 1997, Section H).The interested reader can refer to Appendix C for a discussion on the Fisher-Lange and an implementation of the model using gemact.
For any cell (i, j), with i + j > T , the Fisher-Lange can be used to determine the expected number of future payments ni,j and the future average claim cost mi,j .The original concept behind the CRMR is found in Savelli and Clemente (2009), where the authors assumed that N i,j is Poisson distributed with mean ni,j and Z h;i,j is gamma distributed with mean mi,j .The coefficient of variation of Z h;i,j is ĉzj , the relative variation of the individual payments in development period j, independent from the accident period.The values for ĉzj are estimated using the individual claims data available to the insurer at the evaluation date.The CRMR is extended in Ricotta and Clemente (2016) to take into account of the variability of the severity parameter estimation (Estimation Variance), in addition to the random fluctuations of the underlying process (Process Variance) (Wüthrich & Merz, 2015, p. 28).This is achieved by considering two structure variables (q and ψ), on claim count and average cost, to describe parameter uncertainty on N i,j and Z h;i,j .The parametric assumptions of the model are summarised in Table 4.For i, j = 0, . . ., J , mi,j and nij are obtained from the computation of the Fisher-Lange.

Illustration lossreserve
In this section, we show an example of the CRMR using the lossreserve module.In this respect, we simulated the claims reserving data sets using the individual claim simulator in Avanzi, Taylor, Wang, and Wong (2021) to generate the upper triangles X (T ) , N (T ) , O (T ) and D (T ) .The R (T ) upper triangle has been simulated using the simulator of Avanzi, Taylor, and Wang (2023).The data are simulated using the same assumptions (on an yearly basis) of the Simple, short tail claims environment shown in Al-Mudafer, Avanzi, Taylor, and Wong (2022).
No inflation is assumed to simplify the forecast of the future average costs.Estimating and extrapolating a calendar period effect in claims reserving is a delicate and complex subject, a more detailed discussions can be found in Kuang, Nielsen, andNielsen (2008a, 2008b); Pittarello, Hiabu, and Villegas (2023).Furthermore, in practice, insurers might require a specific knowledge of the environment in which the agents operate to set a value for the inflation (Avanzi et al., 2023).Here, we want to limit the assumptions for this synthetic scenario.4: Parametric assumptions of the CRMR.In the upper table we show the parameters for Z h;i,j , ψ, and q that are gamma distributed.The parameters of the structure variables are specified from the user starting from the variance of ψ and q, indeed Ricotta and Clemente ( 2016) assume E [ψ] = E [q] = 1.The estimator for the average cost of the individual payments is derived with the Fisher-Lange.The variability of the individual payments is instead obtained from the company database.In the lower table we show the parameter for the claim payment number, i.e. a mixed Poisson-gamma distribution with ni,j derived from the Fisher-Lange.

ReserveModel AggregateData
LossReserve Since the entire claim history is available from the simulation and we know the (true) value for the future payments, we can use this information to back-test the performance of our model (Gabrielli & Wüthrich, 2019).In particular, we refer to the actual amount that the insurer will pay in future calendar years as the actual reserve, i.e. the amount that the insurer should set aside to cover exactly the future obligations.The CRMR is compared with the chain-ladder method of Mack (1993), hereafter indicated as CHL, from the R package ChainLadder (Gesmann et al., 2022).These data, together with the data sets from Savelli and Clemente (2014) In more detail, the ReservingModel class arguments are: • The tail parameter (boolean) specifying whether the triangle tail is to be modeled.
• The reserving_method parameter, 'crm' in this case.
• The claims_inflation parameter indicating the vector of claims inflation.In this case, as no claims inflation is present, we simply set it to one in all periods.• The mixing frequency and severity parameters mixing_fq_par and mixing_sev_par.In Ricotta and Clemente (2016) the authors discuss the calibration of the structure variable.Without having a context of a real world example we simply set the structure variables to gamma random variables with mean 1 and standard deviation 0.01, as it is a medium-low risk value in the authors' examples.• The coefficients of variation of the individual claim severity computed for each development period.
In particular, for each development period, the insurer can compute the coefficient of variation from the individual observations and save them in the vector czj.• The vector czj of the coefficients of variation of the individual claim severity, for each development period.The insurer can compute the coefficient of variation from the individual observations, for each development period.
The mean reserve estimate for the CRMR can be extracted from the reserve attribute.In a similarly way to LossModel and LossAggregation, the distribution of the reserve can be accessed from the dist property.
The reserve quantiles can be obtained with the ppf method.Output figures are expressed in millions to simplify reading.

Comparison with the chain-ladder
This section compares the CRMR and the CHL with the actual reserve we know from the simulation.For clarity, the reserve of accident period i is defined as R i = T j=T −i+1 X i,j .Table 5 reports the main results.We can see that the CHL reserve estimate is closer than the CRMR reserve estimate to the actual reserve (the "Reserve" column of each block).In addition to the reserve estimate, we provide the mean squared error of prediction (MSEP) (Wüthrich & Merz, 2015, p. 268) of the reserve estimator for each accident period and for the total reserve.In the aforementioned manuscript, the author introduces the MSEP as a measure of risk under the notion in Wüthrich and Merz (2015, p. 169) conditional on the information available in the upper triangle.The MSEP can be decomposed into two components, the irreducible risk arising from the data generation process (process error, PE) and a risk component arising from the model used to calculate the reserve (model error, ME).We can provide a reference value for the true PE by running multiple simulations from the simulator described in the previous section and calculating the standard error of the reserves by accident period.To obtain the PE, we simulated 100 triangles.Net of the ME component, the results in

Conclusions
This paper introduces gemact, a Python package for non-life actuarial modelling based on the collective risk theory.
Our library expands the reach of actuarial sciences within the growing community of Python programming language.It provides new functionalities and tools for loss modelling, loss aggregation and loss reserving, such as the (a, b, 0) and (a, b, 1) classes of distributions, the AEP algorithm and the collective risk model for claims reserving.Hence, it can be applied in different areas of actuarial sciences, including pricing, reserving and risk management.The package has been designed primarily for the academic environment.Nevertheless, its use as support for insurance business specialists in prototypes modelling, business studies and analyses is not to be excluded.
The structure of our package aims to ensure modifiability, maintainability and scalability, as we thought of gemact as an evolving and growing project in terms of introducing features, integrating functionalities, enhancing methodologies, and expanding its scopes.Possible future enhancements could involve the introduction of new probability distribution families, the implementation of supplementary methodologies for the approximation of quantiles of the sum of random variables, and the addition of costing procedures for exotic and non-traditional reinsurance solutions.

Device specifications
All experiments and analyses were run on a computer with an Intel ® Core™ i7-1065G7 CPU processor with 16GB RAM, running at 1.30GHz, on Windows 11 Home edition.

C Claims reserving with the Fisher-Lange
This section briefly introduces the Fisher-Lange approach.The claims reserve is the sum of the (future) payments, forecast as the product between the predicted future average cost and the predicted number of future payments, for each cell.In formula: R = i+j>J mi,j ni,j . (18) The average claim cost in the lower triangle is forecast as the projection of the inflated average claim cost.
m i,j = m J −j,j i+j h=J +1 (1 + δ h ) , (19)    where δ h represents the claims inflation for calendar period h.
As far the number of claims are concerned, this method assumes that the future number of paid claims is related to the percentage of open claims at the evaluation date and to the claims settlement speed.Indeed, at the evaluation date, the lower triangle is estimated as: where [α j ], j = 0, . . ., J − 1, is the vector of open claims given by α j = E [τ i,j ], and for i = 0, . . ., J − 1 and j = 0, . . ., J − i − 1.It is assumed that α J = 1.The claim settlement speed is then computed for each accident year.The settlement speed for accident period J is where d i represents the number of reported claims for accident period i, with i = 0, . . ., J .The formula is corrected for other accident years following the approach in Savelli and Clemente (2014, p. 141).
Similarly to the CRMR described in Section 5, the results for the Fisher-Lange can be computed with the gemact package.Below, we show an example using the simulated data sets from Section 5.
>>> resmodel = ReservingModel( tail=False, reserving_method='fisher_lange', claims_inflation=claims_inflation) Thereafter, the actual computation of the loss reserve is performed within the LossReserve class: >>> from gemact.lossreserveimport LossReserve >>> lossreserve = LossReserve(data=ad, reservingmodel=resmodel) The LossReserve class comes with a summary view of the estimated reserve per each accident period, in a similar way to the LossModel class, the print_loss_reserve method.In Table 9 we report the CRMR results from  9: Reserves by accident period for the CRMR and the Fisher-Lange.We also report the actual reserve.Amounts are shown in millions.
As expected, being the Fisher-Lange the underlying methodology to the CRMR the results for the claims reserve provided from the two approaches are consistent.
On top of this, insights on the behaviour of Fisher-Lange [α j ] and settlement speed v (i) j , for j = 0, . . ., J − 1, can be inspected with the plot_alpha_fl and plot_ss_fl methods.

Figure 1 :
Figure 1: Illustration of the discretisation methods applied to a gamma(a = 5) severity.The graphs compare the original cdf (blue line) and the discretised (red line) cdf for mass dispersal (top left), upper discretisation (top right), lower discretisation (bottom left), and local moment matching (bottom right) methods.No coverage modifiers are present.

Figure 1
Figure 1 illustrates and graphically contrasts the four different discretisation techniques.

Figure 2 :
Figure 2: Class diagram of the lossmodel module.A rectangle represents a class; an arrow connecting two classes indicates that the target class employs the origin class as an attribute.In this case, a LossModel object entails Frequency, Severity and PolicyStructure class instances.These correspond to the frequency model, the severity model and the policy structure, respectively.The latter, in particular, is in turn specified via one or more Layer objects, which include coverage modifiers of each separate policy component.

Figure 3 :
Figure 3: Class diagram of the lossaggregation module.A rectangle represents a class; an arrow connecting two classes indicates that the target class employs the origin class as an attribute.In this case, a LossAggregation object entails Margins and Copula class instances.

Figure 5 :
Figure 5: Example of events of a non-life insurance claim.

Figure 6 :
Figure 6: Class diagram of the lossreserve module.A rectangle represents a class; an arrow connecting two classes indicates that the target class employs the origin class as an attribute.In this case, a LossReserve object entails ReserveModel and AggregateData class instances.
Iteration 1, the simplex S1 in the first iteration is approximated with the square Q1.

Figure 7 :
Figure 7: Sketch of the first two iterations of the AEP algorithm in the two-dimensional case.
The aggregate loss mean, standard deviation and skewness can be accessed with the mean, std, and skewness methods, respectively.The code below shows how to use these methods.
More precisely, in the Layer class the contract cover and the deductible are provided.Once the model assumptions are set and the policy structure is specified, the aggregate loss distribution can be computed.
Coverage modifiers are set in a PolicyStructure object.
Contrarily, once the aggregate loss distribution is determined (via dist_calculate method), all layer premiums in pure_premium_dist are available from the costing (costing method).As expected, pure_premium content remains unaffected.

Table 2 :
Comparison for different contracts of aggregate and gemact implementation of the aggregate loss distribution computation via FFT.When there are no aggregate coverage modifiers, reference values are given in addition to estimated ones.For the XL and the XL w/agg.contracts, individual conditions are c = 1000000 and d = 10000; for the SL and the XL w/agg.contracts, aggregate coverage modifiers are u = 1000000 and v = 50000.Execution times are expressed in seconds.
, have been saved in the gemdata module for reproducibility.The lossreserve class diagram is illustrated in Figure6: the triangular data sets are stored in the AggregateData class, and the model parameters are contained in the ReservingModel class.The actual computation of the loss reserve is performed with the LossReserve class that takes as inputs an AggregateData object, a ReservingModel object and the parameters of the claims reserving model computation.

Table 5 :
Table5appear to show that the CHL underestimates the forecast error in the early accident periods and overestimates the forecast error in the later accident periods.Reserve and mean squared error of prediction (MSEP) by accident period for the CRMR and the CHL.The actual reserve and its process error (PE) by accident period are also indicated.Amounts are shown in millions.To conclude, Table6shows the total reserves of the CRMR and the CHL.In particular, it highlights the reserve amount as a percentage of the actual reserve (8599.04)and the coefficient of variation (CoV), i.e. the MSEP divided by the reserve.The total relative variability referred to the PE component gives a CoV of 0.76 %.

Table 6 :
Total reserve estimates, their relative value, as a fraction of the actual value (8599.04),and their coefficients of variation (CoV), for the CRMR and the CHL.Absolute amounts are reported in millions.
Table 5 and we add the results for the Fisher-Lange.