3.1 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, Reference Parodi2014, 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:

  (4) \begin{equation} \text{P}^{QS} = \mathbb{E}\left [ Q_a \left ( X \right)\right]. \end{equation}

• • 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 \geq 0$ :

  (5) \begin{equation} \text{P}^{XL} = \mathbb{E}\left [ \sum _{i=1}^{N} L_{c,d} (Z_i) \right]. \end{equation}

• • 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 \geq 0$ :

  (6) \begin{equation} \text{P}^{SL} = \mathbb{E}\left [ L_{u, v} (X) \right]. \end{equation}

The model introduced by Equation (1) and implemented in gemact can be used for costing contracts like the XL with Reinstatements (RS) in Sundt (Reference Sundt1990). Assuming the aggregate cover $u$ is equal to $ (K + 1) c$ , with $K \in \mathbb{Z}^+$ :

(7) \begin{equation} \text{P}^{RS} = \frac{\mathbb{E}\left [ L_{u, v} (X) \right]}{1+\frac{1}{c} \sum _{k=1}^K l_k \mathbb{E}\left [ L_{c, (k-1)c+v}(X) \right]}, \end{equation}

where $K$ is the number of reinstatement layers and $l_k \in [0, 1]$ is the reinstatement premium percentage, with $k=1, \ldots, 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, that is, $l_k \text{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, Reference Parodi2014, 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 $\text{P}^{RS}$ is beyond the scope of this manuscript; the interested reader can refer to Sundt (Reference Sundt1990), Parodi (Reference Parodi2014, p. 325), and Antal (Reference Antal2009, p. 57).

3.2 Computational methods for the aggregate loss distribution

The cumulative distribution function (cdf) of the aggregate loss in Equation (1) is

(8) \begin{align} F_X(x)=P[X \leq x]=\sum _{k=0}^{\infty } p_k F_Z^{* k}(x) \end{align}

where $p_k=P[N=k]$ , $F_Z(x)=P\left [ Z \leq x \right]$ and $F_Z^{* k}(x)=P\left [ Z_1 + \ldots + Z_k \leq x \right]$ .

Moreover, the characteristic function of the aggregate loss $\phi _{X}\;:\; \mathbb{R} \rightarrow \mathbb{C}$ can be expressed in the form:

(9) \begin{equation} \phi _{X}(t)=\mathcal{P}_{N}\left (\phi _{Z}(t)\right), \end{equation}

where $\mathcal{P}_{N}(t)=\mathrm{E}\left [t^N\right]$ is the probability generating function of the frequency $N$ and $\phi _{Z}(t)$ is the characteristic function of the severity $Z$ (Klugman et al., Reference Klugman, Panjer and Willmot2012, p. 153).

The distribution in Equation (8), except in a few cases, cannot be computed analytically, and its direct calculation is numerically expensive (Parodi, Reference Parodi2014, p. 239). For this reason, different approaches have been analyzed to approximate the distribution function of the aggregate loss, including parametric and numerical quasi-exact methods (for a detailed treatment refer to Shevchenko, Reference Shevchenko2010). Among the latter, gemact implements MC simulation (Klugman et al., Reference Klugman, Panjer and Willmot2012, p. 467), DFT (Bühlmann, Reference Bühlmann1984; Grübel & Hermesmeier, Reference Grübel and Hermesmeier1999; Wang, 1998), and the so-called recursive formula (Panjer, Reference Panjer1981). A brief comparison of accuracy, flexibility, and speed of these methods can be found in Parodi (Reference Parodi2014, p. 260) and Wüthrich (Reference Wüthrich2023, p. 127). This section details these last two computational methods based on discrete mathematics to approximate the aggregate loss distribution.

Henceforth, let us consider, for $j=0,1,2,\ldots,m-1$ and $h\gt 0$ , an arithmetic severity distribution with probability sequence:

\begin{align*} \textbf{\{f\}}=\{f_0, f_1, \ldots, f_{m-1}\}, \end{align*}

where $f_j=P[Z=j\cdot h]$ . The discrete version of Equation (8) becomes

\begin{equation*} g_s=\sum _{k=0}^{\infty } p_k f_s^{* k}, \end{equation*}

where $g_s = P[X=s]$ and

\begin{equation*} f_s^{* j}\;:\!=\; \begin {cases}1 & \text { if } k=0 \text { and } s=0 \\[5pt] 0 & \text { if } k=0 \text { and } s \in \mathbb {N} \\[5pt] \sum _{i=0}^s f_{s-i}^{*(k-1)} f_i & \text { if } k \gt 0.\end {cases} \end{equation*}

3.2.1 Discrete fourier transform

The DFT of the severity $\textbf{\{f\}}$ is, for $k=0,\ldots,m-1$ , the sequence

\begin{equation*} \{\hat{\textbf{f}}\}= \{\hat{f}_0, \hat {f}_1, \ldots, \hat{f}_{m-1} \},\end{equation*}

where

(10) \begin{equation} \widehat{f}_k=\sum _{j=0}^{m-1}f_j e^{\frac{2\pi i k j}{m}}. \end{equation}

The original sequence can be reconstructed with the inverse DFT:

\begin{equation*}f_j=\frac {1}{m} \sum _{k=0}^{m-1} \widehat {f_k} e^{-\frac {2 \pi i k j}{ m}}.\end{equation*}

The sequence of probabilities $\textbf{\{g\}}=\{g_0, g_1, \ldots, g_{m-1}\}$ can be approximated taking the inverse DFT of

(11) \begin{equation} \{ \hat{\textbf{g}}\}\;:\!=\; \mathcal{P}_{N}\left (\{\hat{\textbf{f}}\}\right). \end{equation}

The original sequence can be computed efficiently with a FFT algorithm, when $m$ is a power of $2$ (Embrechts & Frei, Reference Embrechts and Frei2009).

3.2.2 Recursive formula

Assume that the frequency distribution belongs to the $(a, b, 0)$ class, that is, for $k \geq 1$ and $a, b\in \mathbb{R}$ :

(12) \begin{equation} p_k=\left (a+\frac{b}{k}\right)p_{k-1}. \end{equation}

Here, $p_0$ is an additional parameter of the distribution (Klugman et al., Reference Klugman, Panjer and Willmot2012, p. 505). The $(a, b, 0)$ class can be generalized to the $(a, b, 1)$ class assuming that the recursion in Equation (12) holds for $k=2,3,\ldots$

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 $\textbf{{g}}$ can obtained recursively using the following formula:

(13) \begin{equation} g_s=\frac{\left [p_1-(a+b) p_0\right] f_s+\sum _{j=1}^s(a+b j/ s) f_j g_{s-j}}{1-a f_0}, \end{equation}

with $1 \leq s \leq m-1$ and given the initial condition $g_0=\mathcal{P}_{N}\left (f_0\right)$ .

3.3 Severity discretization

The calculation of the aggregate loss with DFT or with the recursive formula requires an arithmetic severity distribution (Embrechts & Frei, Reference Embrechts and Frei2009). Conversely, the severity distribution in Equation (1) is often calibrated on a continuous support. In general, one needs to choose a discretization approach to arithmetize the original distribution. This section illustrates the methods for the discretization of a continuous severity distribution available in the gemact package.

Let $F_z\;:\; \mathbb{R^+} \rightarrow [0, 1]$ be the cdf of the distribution to be discretized. For a bandwidth, or discretization step, $h\gt 0$ and an integer $m$ , a probability $f_j$ is assigned to each point $h j$ , with $j= 0, \ldots,m-1$ . Four alternative methods for determining the values for $f_j$ are implemented in gemact.

1. 1. The method of mass dispersal:

   \begin{align*} f_j= \begin{cases}F_Z \left ( \frac{h}{2} \right) & j=0 \\[5pt] F_Z\left ( hj+ \frac{h}{2} \right)-F_Z\left (hj-\frac{h}{2}\right) & j=1,\ldots,m-2\\[5pt] 1-F_Z\left (hj-\frac{h}{2}\right) & j=m-1. \end{cases} \end{align*}

2. 2. The method of the upper discretization:

   \begin{align*} f_j= \begin{cases}F_Z\left ( hj+ h \right)-F_Z\left (hj\right) & j=0,\ldots,m-2\\[5pt] 1-F_Z\left (hj\right) & j=m-1. \end{cases} \end{align*}

3. 3. The method of the lower discretization:

   \begin{align*} f_j= \begin{cases}F_Z \left ( 0 \right) & j=0 \\[5pt] F_Z\left ( hj\right)-F_Z\left (hj-h\right) & j=1,\ldots,m-1.\\[5pt] \end{cases} \end{align*}

4. 4. The method of local moment matching:

   \begin{align*} f_j= \begin{cases}1-\frac{\mathbb{E}[Z\wedge h] }{h} & j=0 \\[5pt] \frac{2 \mathbb{E}[Z \wedge h j]-\mathbb{E}[Z \wedge h (j-1)]-\mathbb{E}[Z \wedge h (j+1)]}{h} & j=1,\ldots,m-1\end{cases} \end{align*}
   where $\mathbb{E}[Z\wedge h] = \int _{-\infty }^h t \; dF_z(t)+h[1-F_z(h)]$ .

Fig. 1 illustrates and graphically contrasts the four different discretization techniques.

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

The default discretization method in gemact is the method of the mass dispersal, in which each $h j$ point is assigned with the probability mass of the $h$ -span interval containing it (Fig. 1 top left). The upper discretization and lower discretization 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 Fig. 1, top right and bottom left graphs). The method of local moment matching (Fig. 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 (Reference Gerber1982). 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, Reference Embrechts and Frei2009). 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 discretization 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 discretization 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 discretization parameters can be found in Parodi (Reference Parodi2014, 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, Reference Embrechts and Frei2009). The reader should refer to Klugman et al. (Reference Klugman, Panjer and Willmot2012, p. 179) and Embrechts & Frei (Reference Embrechts and Frei2009) for a more detailed treatment and discussion of discretization methods.

The above-mentioned discretization methods are modified accordingly to reflect the cases where the transformation of Equation (3) is applied to the severity (Klugman et al., Reference Klugman, Panjer and Willmot2012, p. 517). Below, we illustrate how to perform severity discretization in gemact.

Once a continuous distribution is selected from those supported in our package (see Section A), the severity distribution is defined via the Severity class.

>>> from gemact.lossmodel import 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.

>>> massdispersal = severity.discretize(

discr_method='massdispersal',

n_discr_nodes = 50000,

discr_step = .01,

deductible = 0

)

In order to perform the discretization, the following arguments are needed:

• • The chosen discretization method via discr_method.

• • The number of nodes ( $m$ ) set in the n_discr_nodes argument.

• • The severity discretization step ( $h$ ) is in the discr_step argument.

• • If necessary, a deductible specifying where the discretization begins. The default value is zero.

After the discretization is achieved, the mean of the discretized 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 discretization can be visually examined using the plot_discretized_severity_cdf. This method is based on the pyplot interface to matplotlib (Hunter, Reference Hunter2007). 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.

>>> from gemact import distributions

>>> dist = distributions.Gamma(a = 5)

>>> nodes = np.arange(0, 20)

>>> true_cdf = dist.cdf(nodes)

Next, we plot the discrete severity using the plot_discr_sev_cdf method.

>>> import matplotlib.pyplot as plt

>>> severity.plot_discr_sev_cdf(

discr_method='massdispersal',

n_discr_nodes = 20,

discr_step = 1,

deductible = 0,

color='#a71429'

)

>>> plt.plot(nodes, true_cdf, color='#4169E1')

>>> plt.title('Method of Mass Dispersal')

>>> plt.xlabel('z')

>>> plt.show()

The arguments discr_method, n_discr_nodes, discr_step, and deductible can be used in the same manner as those described in the discretize method. The argument color is from the matplotlib.pyplot.plot method. The methods title and xlabel were also exported from matplotlib.pyplot.plot to add custom labels for the title and the x-axis (Hunter, Reference Hunter2007).

The output of the previous code is shown in the top left graph of Fig. 1. To obtain the other graphs, simply set discr_method to the desired approach, that is, ‘upper_discretisation,’ ‘lower_discretisation’, or ‘localmoments.’

3.4 Supported distributions

The gemact package makes for the first time the $(a, b, 0)$ and $(a, b, 1)$ distribution classes (Klugman et al., Reference Klugman, Panjer and Willmot2012, 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 MC simulation, with the random generator method for the ZTPoisson class.

>>> random_variates = ztpois.rvs(10**5, random_state = 1)

>>> print('Simulated Mean: ', np.mean(random_variates))

Simulated Mean: 2.3095

>>> print('Exact Mean: ', ztpois.mean())

Exact Mean: 2.3130352854993315

Furthermore, supported copula functions can be accessed via the copulas module. Below, we compute the cdf of a two-dimensional Gumbel copula.

>>> 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 and 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 Section 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 modeling in general. An example of this is the discretize method. A similar reasoning applies to the Frequency class.

3.5 Illustration lossmodel

The following are examples of how to get started and use lossmodel module and its classes for costing purposes. As an overview, Fig. 2 schematizes 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.

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.

>>> from gemact.lossmodel import 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, Reference Parodi2014, 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 Subsection 3.3.

>>> from gemact.lossmodel import LossModel

>>> lm_mc = LossModel(

frequency=frequency,

severity=severity,

aggr_loss_dist_method='mc',

n_sim = 10**5,

random_state = 1

)

INFO:lossmodel|Approximating aggregate loss distribution via Monte Carlo simulation

INFO:lossmodel|MC simulation completed

In the previous example, in more detail, lm_mc object adopts the MC 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 as follows:

• • the number of simulations n_sim,

• • the (pseudo)random number generator initializer 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.

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.

>>> lm_mc.mean(use_dist=True)

19.963155575580227

>>> lm_mc.mean(use_dist=False)

20.0

>>> lm_mc.coeff_variation(use_dist=True)

0.5496773171375182

>>> lm_mc.coeff_variation(use_dist=False)

0.5477225575051661

>>> lm_mc.skewness(use_dist=True)

0.6410913579225725

>>> lm_mc.skewness(use_dist=False)

0.6390096504226938

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 closed-form 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, Reference Parodi2014, 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. Coverage modifiers are set in a PolicyStructure object.

>>>from gemact.lossmodel import PolicyStructure, Layer

>>> policystructure = PolicyStructure(

layers=Layer(

cover = 20,

deductible = 5

))

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.

>>> lm_XL = LossModel(

frequency=frequency,

severity=severity,

policystructure=policystructure,

aggr_loss_dist_method='fft',

sev_discr_method='massdispersal',

n_aggr_dist_nodes = 2**17

)

INFO:lossmodel|Approximating aggregate loss distribution via FFT

INFO:lossmodel|FFT completed

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 discretize the severity distribution sev_discr_method. Above, we opted for the method of mass dispersal ("massdispersal").

• • The number of nodes of the discretized severity n_sev_discr_nodes (optional).

• • The discretization step sev_discr_step (optional). When a cover is present, gemact automatically adjusts the discretization step parameter to have the correct number of nodes in the transformed severity support.

The same arguments shall be specified while computing the aggregate loss distribution with the recursive formula, that is, aggr_loss_dist_method set to "recursive."

The costing specifications of a LossModel object can be accessed with the method print_costing_specs().

The previous output exhibits a summary of the contract structure (cover, deductible, aggregate cover, and aggregate deductible) and details about the costing results. It is noted that the default value of the share participation equals 1, that is, $a$ is equal to $1$ in Equation (2).

Similar to the previous example, user can access moments of the aggregate loss calculated from the approximated distribution and from the closed-from solution.

>>> lm_XL.mean(use_dist=True)

3.509346100359707

>>> lm_XL.mean(use_dist=False)

3.50934614394912

>>> lm_XL.coeff_variation(use_dist=True)

1.0001481880266856

>>> lm_XL.coeff_variation(use_dist=False)

1.0001481667319252

>>> lm_XL.skewness(use_dist=True)

1.3814094240544392

>>> lm_XL.skewness(use_dist=False)

1.3814094309741256

The next example illustrates the costing of an XL with RS. The PolicyStructure object is set as follows.

>>> policystructure_RS = PolicyStructure(

layers=Layer(

cover = 100,

deductible = 0,

aggr_deductible = 100,

reinst_percentage = 1,

n_reinst = 2

))

The relevant parameters are as follows:

• • the aggregate deductible parameter aggr_deductible,

• • the number of reinstatements n_reinst,

• • the reinstatement percentage reinst_percentage.

Below, we compute the pure premium of Equation (7), given the parametric assumptions on the frequency and the severity of the loss model.

>>> lm_RS = LossModel(

frequency=Frequency(

dist='poisson',

par={'mu': .5}

),

severity=Severity(

dist='pareto2',

par={'scale': 100, 'shape': 1.2}

),

policystructure = policystructure_RS,

aggr_loss_dist_method='fft',

sev_discr_method='massdispersal',

n_aggr_dist_nodes = 2**17

)

>>> print(Pure premium (RS):, lm_RS.pure_premium_dist[0])

Pure premium (RS): 4.319350355177216

A PolicyStructure object can handle multiple independent 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, and the third with an aggregate cover.

>>> policystructure=PolicyStructure(

layers=[

Layer(cover = 100, deductible = 100, share = 0.5),

Layer(cover = 200, deductible = 100, n_reinst = 2,

reinst_percentage = 0.6),

Layer(cover = 100, deductible = 100, aggr_cover = 200)

])

>>> lossmodel_multiple = LossModel(

frequency=Frequency(

dist='poisson',

par={'mu': .5}

),

severity=Severity(

dist='genpareto',

par={'loc': 0, 'scale': 83.34, 'c': 0.834}

),

policystructure=policystructure

)

WARNING:lossmodel|Aggregate loss distribution calculation is omitted as

aggr_loss_dist_method is missing

WARNING:lossmodel|Layer 2: costing is omitted as aggr_loss_dist_method is missing

WARNING:lossmodel|Layer 3: costing is omitted as aggr_loss_dist_method is missing

As outlined by the warning messages, since the instantiation of lossmodel_multiple lacks of aggr_loss_dist_method, the calculation of the aggregate loss distribution is omitted. Therefore, costing results are accessible solely for the first Layer, since this is the only one without aggregate coverage modifiers. This fact is reflected in pure_premium and pure_premium_dist properties, containing the premiums derived from the closed-form means and the approximated aggregate loss distribution means, respectively. The latter are indeed not available.

>>> lossmodel_multiple.pure_premium

[8.479087307840043, None, None]

>>> lossmodel_multiple.pure_premium_dist

[None, None, None]

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.

>>> lossmodel_multiple.dist_calculate(

aggr_loss_dist_method='fft',

sev_discr_method='massdispersal',

n_aggr_dist_nodes = 2**17

)

INFO:lossmodel|Computation of layers started

INFO:lossmodel|Computing layer: 1

INFO:lossmodel|Approximating aggregate loss distribution via FFT

INFO:lossmodel|FFT completed

INFO:lossmodel|Computing layer: 2

INFO:lossmodel|Approximating aggregate loss distribution via FFT

INFO:lossmodel|FFT completed

INFO:lossmodel|Computing layer: 3

INFO:lossmodel|Approximating aggregate loss distribution via FFT

INFO:lossmodel|FFT completed

INFO:lossmodel|Computation of layers completed

>>> lossmodel_multiple.costing()

>>> lossmodel_multiple.pure_premium_dist

[8.479087307062226, 25.99131088702302, 16.88704720494799]

>>>lossmodel_multiple.pure_premium

[8.479087307840043, None, None]

At last, it is remarked that each Layer in PolicyStructure is associated with an index idx, starting from 0, based on the layer order of the instantiation of the PolicyStructure object. This is of help when the user needs to retrieve particular information and features, or to apply methods to one specific layer. All the methods that include the idx argument have 0 as default value, meaning that they are applied to the first (or only) layer unless otherwise specified. For example, the print_policy_layer_specs method produces a table of recap of the features of the layer indicated by the idx argument. Below idx equals 1, namely the second Layer in policystructure.

Likewise, the code block below returns the aggregate loss distribution mean of the third Layer, as idx is set to 2.

>>> lossmodel_multiple.mean(idx = 2)

16.88704720494799

3.6 Comparison of the methods for computing the aggregate loss distribution

In this section, we analyze accuracy and speed of the computation of the aggregate loss distribution using FFT, recursive formula (recursion), and MC simulation approaches, as the number of nodes, the discretization 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 (Reference Parodi2014, 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.68$ and $177728.30$ , respectively.

• • Frequency: Poisson distribution with parameter $\mu = 3$ , with analysis threshold $d$ . It belongs to the $(a, b, 0)$ family described in Subsection 3.2.2 with parameters $a=0$ , $b=3$ and $p_0=e^{-3}$ .

• • Policy structure: contract with deductible $d=10000$ .

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, Reference Bean2000, p. 382):

\begin{align*} & \mathbb{E}\left [X\right] = \mu \mathbb{E}\left [ L_{d, \infty } (Z) \right],\\[5pt] & \text{CoV}\left [X\right] = \dfrac{ \mathbb{E}\left [ L_{d, \infty } (Z)^2 \right]^{1/2}}{\mu ^{1/2} \mathbb{E}\left [ L_{d, \infty } (Z) \right]},\\[5pt] & \text{Skewness}\left [X\right] = \dfrac{ \mathbb{E}\left [ L_{d, \infty } (Z)^3 \right]}{\mu ^{1/2} \mathbb{E}\left [ L_{d, \infty } (Z)^2 \right]^{3/2}}. \end{align*}

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 et al., Reference Martelli, Ravenscroft and Ascher2005, 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. We considered for FFT and recursion the values 50, 100, 200, and 400 for the discretization step $h$ , and $2^{14}, 2^{16}, 2^{18}$ for the number of nodes. It can be noted that FFT and recursion produce similar figures in terms of accuracy, but the former is drastically faster. This is expected as FFT takes essentially $\mathcal{O}(m \log (m))$ operations, compared to the $\mathcal{O}(m^2)$ operations for recursion (Embrechts & Frei, Reference Embrechts and Frei2009). 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.

3.7 Comparison with aggregate FFT implementation

The aggregate package (Mildenhall, Reference Mildenhall2022) allows to compute the aggregate loss distribution using FFT. In this section, we compare our implementation with aggregate (version 0.9.3) implementation to show that the two provide similar results. We adopt the same underlying frequency and severity assumptions of Subsection 3.6 and contracts with different combinations of individual and aggregate coverage modifiers. In particular, we first consider no reinsurance, then an XL, with individual-only coverage modifiers, a SL, and finally an XL with individual and aggregate coverage modifiers (XL w/agg.). In line with the example in Parodi (Reference Parodi2014, p. 262), individual coverage modifiers are $c= 1 000 000$ and $d= 10 000$ , and aggregate coverage modifiers are $u= 1 000 000$ and $v= 50 000$ . The number of nodes $m$ is set to $2^{22}$ in all the calculations.

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 (MC) simulation when varying the number of nodes ( $m$ ), the discretization 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

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 initialization 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., Reference Martelli, Ravenscroft and Ascher2005, 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.

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= 1 000 000$ and $d= 10 000$ ; for the SL and the XL w/agg. contracts, aggregate coverage modifiers are $u= 1 000 000$ and $v= 50 000$ . Execution times are expressed in seconds

4.1 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 Fig. 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, that is, the copula.

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.

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 A.1 and Table A.2 of Section A for the complete list of the supported distributions and copulas.

>>> from gemact import LossAggregation, Copula, Margins

>>> lossaggregation = LossAggregation(

margins=Margins(

dist=['genpareto', 'lognormal'],

par=[{'loc': 0, 'scale': 1/.9, 'c': 1/.9}, {'loc': 0, 'scale': 10,

'shape': 1.5}],

),

copula=Copula(

dist='frank',

par={'par': 1.2, 'dim': 2}

),

n_sim = 500000,

random_state = 10,

n_iter = 8

)

Besides marginal and copula assumptions, the instantiation of LossAggregation accepts the arguments n_sim and random_state to set the number of simulation and the (pseudo)random number generator initializer of the MC simulation. If n_sim and the random_state are omitted, the execution is bypassed. Nonetheless, the user can perform it at a later point through the dist_calculate method. Also, n_iter controls the number of iterations of the AEP algorithm. This parameter is optional (default value is 7) and can be specified either at the creation of class or directly in methods that include the use of this approach.

In the following code block, we show how to calculate the cdf of the sum of a generalized Pareto and a lognormal-dependent random variables, using the AEP algorithm ("aep") and the MC simulation approach ("mc"). The underlying dependency structure is a Frank copula.

>>> s = 300 # arbitrary value

>>> p_aep = lossaggregation.cdf(x = s, method=aep)

>>> print('P(X1+X2 <= s) = ', p_aep)

P(X1+X2 <= s) = 0.9811620158197308

>>> p_mc = lossaggregation.cdf(x = s, method='mc')

>>> print('P(X1+X2 <= s) = ', p_mc)

P(X1+X2 <= s) = 0.98126

LossAggregation includes other functionalities like the survival function sf, the quantile function ppf, and the generator of random variates rvs. Furthermore, for the MC simulation approach, it is possible to derive empirical statistics and moments using methods such as for example moment, mean, var, skewness, and censored_moment. To conclude, the last code block illustrates the calculation of the quantile function via the ppf method.

>>> lossaggregation.ppf(q = p_aep, method=aep)

300.0000003207744

>>> lossaggregation.ppf(q = p_mc, method=mc)

299.9929982860278

4.2 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 MC 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. (Reference Arbenz, Embrechts and Puccetti2011) 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 $\gamma$ of the marginal distributions are $0.9$ and $1.8$ ; the Clayton copula has parameter $\theta =1.2$ . In three dimensions ( $d=3$ ), the additional marginal has parameter $\gamma =2.6$ , and the copula has $\theta =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., Reference Martelli, Ravenscroft and Ascher2005, 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.

Table 3. Accuracy and speed of the cdf calculation using the AEP algorithm and the Monte Carlo (MC) simulation approach 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. (Reference Arbenz, Embrechts and Puccetti2011). 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 (labeled with a ^*) lists those of the original manuscript

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 MC 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 MC 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 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. Fig. 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 $\rho$ (the non-diagonal entry of the correlation matrix) and $\theta$ parameters, for the Gaussian copula and the Clayton copula, respectively. The cdf values were also compared with those obtained by the MC 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.

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, and 1.95 for the bivariate Gaussian copula (left plot) and 2, 2.3, 2.65, and 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 $\rho$ and $\theta$ change. The results of the AEP algorithm correspond to those of the Monte Carlo (MC) simulation approach, 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.

5.1 Problem framework

Let the index $i=0,\ldots,\mathcal{J}$ denotes the claim accident period, and let the index $j=0,\ldots,\mathcal{J}$ represents the claim development period, over the time horizon $\mathcal{J}\gt 0$ . The so-called development triangle is the set:

\begin{equation*} \mathcal {T}=\{(i, j)\;: \; i=0, \ldots, \mathcal {J}, \; j=0, \ldots, \mathcal {J}; \; i+j \leq \mathcal {J}\}, \quad \mathcal {J}\gt 0. \end{equation*}

Below is the list of the development triangles used in this section.

• • The triangle of incremental paid claims:

  \begin{equation*}X^{(\mathcal {T})} =\left \{x_{i, j}\;:\; (i, j) \in \mathcal {T} \right \}, \end{equation*}
  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:

  \begin{equation*}N^{(\mathcal {T})} =\left \{n_{i, j}\;:\; (i, j) \in \mathcal {T} \right \}, \end{equation*}
  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
  \begin{equation*}m_{i,j}=\frac {x_{i,j}}{n_{i, j}},\end{equation*}
  where $m_{i,j}$ is the average claim cost for accident period $i$ and development period $j$ .

• • The triangle of incremental amounts at reserve:

  \begin{equation*}R^{(\mathcal {T})} =\left \{r_{i, j}\;:\; (i, j) \in \mathcal {T} \right \}, \end{equation*}
  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:

  \begin{equation*}O^{(\mathcal {T})} =\left \{o_{i, j}\;:\; (i, j) \in \mathcal {T} \right \},\end{equation*}
  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:

  \begin{equation*}D^{(\mathcal {T})} =\left \{d_{i, j}\;:\; (i, j) \in \mathcal {T} \right \}, \end{equation*}
  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 & Clemente (Reference Ricotta and Clemente2016) and Clemente et al. (Reference Clemente, Savelli and Zappa2019), hereafter referred to as CRMR, which connects claims reserving with aggregate distributions. Indeed, the claims reserve is the sum of the (future) payments:

(16) \begin{equation} R = \sum _{i+j\gt \mathcal{J}}X_{i,j}. \end{equation}

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 (Reference Ricotta and Clemente2016). For $i,j=0,\ldots,\mathcal{J},$

(17) \begin{equation} X_{i,j}=\sum _{h=1}^{N_{i,j}} \psi Z_{h;i,j}, \end{equation}

where the index $h$ is referred to the individual claim severity. The random variable $\psi$ follows a gamma distribution, see Table 4, and introduces dependence between claim-sizes of different cells.

Table 4. Parametric assumptions of the CRMR. In the upper table, we show the parameters for $Z_{h;i,j}, \psi,$ and $q$ that are gamma-distributed. The parameters of the structure variables are specified from the user starting from the variance of $\psi$ and $q$ , indeed Ricotta and Clemente (Reference Ricotta and Clemente2016) assume $\mathbb{E}[\psi]=\mathbb{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, that is, a mixed Poisson-gamma distribution with $\hat n_{i,j}$ derived from the Fisher–Lange

5.2 Illustration lossreserve

In this section, we show an example of the CRMR using the lossreserve module. In this respect, we simulated the claims reserving datasets using the individual claim simulator in Avanzi et al. (Reference Avanzi, Taylor, Wang and Wong2021) to generate the upper triangles $X^{\left (\mathcal{T}\right)}, N^{\left (\mathcal{T}\right)}, O^{\left (\mathcal{T}\right)}$ and $D^{\left (\mathcal{T}\right)}$ . The $R^{\left (\mathcal{T}\right)}$ upper triangle has been simulated using the simulator of Avanzi et al. (Reference Avanzi, Taylor and Wang2023). The data are simulated using the same assumptions (on an yearly basis) of the Simple, short tail claims environment shown in Al-Mudafer et al. (Reference Al-Mudafer, Avanzi, Taylor and Wong2022).

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, and a more detailed discussions can be found in Kuang et al. (Reference Kuang, Nielsen and Nielsen2008a, Reference Kuang, Nielsen and Nielsenb), Pittarello et al. (Reference Pittarello, Hiabu and Villegas2023). 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., Reference Avanzi, Taylor and Wang2023). Here, we want to limit the assumptions for this synthetic scenario.

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, Reference Gabrielli and Wüthrich2019). In particular, we refer to the actual amount that the insurer will pay in future calendar years as the actual reserve, that is, 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 (Reference Mack1993), hereafter indicated as CHL, from the R package ChainLadder (Gesmann et al., Reference Gesmann, Murphy, Zhang, Carrato, Wüthrich, Concina and Dal Moro2022). These data, together with the datasets from Savelli & Clemente (Reference Savelli and Clemente2014), have been saved in the gemdata module for reproducibility.

>>> from gemact import gemdata

>>> ip = gemdata.incremental_payments_sim

>>> pnb = gemdata.payments_number_sim

>>> cp = gemdata.cased_payments_sim

>>> opn = gemdata.open_number_sim

>>> reported = gemdata.reported_claims_sim

>>> czj = gemdata.czj_sim

The lossreserve class diagram is illustrated in Fig. 6: the triangular datasets 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.

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.

>>> from gemact import AggregateData

>>> from gemact import ReservingModel

>>> import numpy as np

>>> ad = AggregateData(

incremental_payments=ip,

cased_payments=cp,

open_claims_number=opn,

reported_claims=reported,

payments_number=pnb

)

>>> resmodel_crm = ReservingModel(

tail=False,

reserving_method='crm',

claims_inflation=np.array([1]),

mixing_fq_par = .01,

mixing_sev_par = .01,

czj=czj

)

In more detail, the ReservingModel class arguments are as follows:

• • 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 & Clemente (Reference Ricotta and Clemente2016), 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 CoV 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 CoV from the individual observations, for each development period.

The computation of the loss reserve occurs within the LossReserve class:

>>> from gemact import LossReserve

>>> lr = LossReserve(

data=ad,

reservingmodel=resmodel_crm,

ntr_sim = 1000,

random_state = 42

)

To instantiate the LossReserve class, AggregateData and ReservingModel objects need to be passed as arguments. The additional arguments required for the computation are as follows:

• • the number of simulated triangles ntr_sim,

• • the (pseudo)random number generator initializer random_state.

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.

>>> lr.reserve

8245996481.515498

>>> lr.ppf(q = np.array([.25, .5, .75, .995, .9995]))/10**6

array([8156.79994049, 8244.66099324, 8335.94438221, 8600.42263791, 8676.84206101])

Table 5. 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

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