Stochastic Claims Reserve in the Healthcare System: A New Methodology Applied in Italy

1. Introduction

1.1. The medical malpractice in Italy: liability rule and definition of claims reserve

The term “medical malpractice” is intended to define those damages (tending to be biological, but also in perspective moral) arising from errors in the health sector of the practitioner and closely related to clinical activity in the strict sense [1]. This damage is closely correlated with civil liability in healthcare, which profoundly impacts the economy and society. Indeed, the healthcare expenditure was estimated for 2019 at 8.7% of the national Gross Domestic Product (GDP), which is below the EU average of 9.9% [3]. It should be noted that approximately 15% of hospital expenditure in high-income countries is used to correct complications of care and patient harm [2]. Although in these countries with an efficient healthcare system, adverse events seem to be progressively decreasing. It is estimated that 10% of patients experience an adverse event during medical-surgical treatment, resulting in a malpractice injury. The most frequently reported types of adverse events, on a global level, undoubtedly include surgical procedures - with significant peaks among orthopaedics and traumatology, general surgery and obstetrics and gynaecology [6] - diagnostic errors, accidental falls, and the contraction of nosocomial infections, in reference to which specific investigation would be due. In Italy, no official national published statistics on medical adverse events and/or claims exist. However, some regions have published independent reports on hospital claims in the past, showing that most complaints are received in the above-mentioned operational areas [18,24]. When an adverse event occurs, Italian healthcare liability is imposed on healthcare organizations and, potentially, healthcare professionals who committed the error. Healthcare public and private institutions are liable with contractual liability for the behaviour of healthcare professionals whose services they employ, with reference to professionals, including researchers, under any contract of the National Health Service. For this reason, public and private health and social care facilities must be provided with insurance coverage or other similar measures for third parties’ liability and towards employees. Other measures include the possibility for organizations to constitute cover funds from internal provisions to protect all or part of the health liability risk. The amount of money to pay insurance premiums or to provide coverage funds derives from the national health fund, which is then distributed to the different regions. In turn, they may distribute the money to their health organizations or centralize the resources to shield against claims. Therefore, insurance or self-insurance coverage must be provided to be operative in the event of a medical malpractice claim. However, given the significant fluctuation in medical malpractice claims, the challenge arises regarding how much money should be allocated annually to guarantee the fulfilment of its financial commitments and the payment of compensation for damages. In this study, the authors aim to address the following research questions: 1) Which models are appropriate for estimating claims reserves in the healthcare sector? 2) How can the models be adapted to enhance the state-of-the-art (SoA) models to fit our specific context?

The proposed work can be positioned in the context described above to build a solid procedure to estimate and manage the claims reserve. The treatment developed here will be technical because the authors want to provide a complete mathematical and computational tool for the healthcare system to allocate appropriate resources. The paper is organized as follows: 1) the first part is focused on the mathematical background behind the GLM to underline the foundation of the applied methodology, jointly with the definition of the quasi-likelihood function to model the response variables, as well as the algorithms to locate the best parameters that maximize these function and the corresponding errors associated to the parameters of the quasi-likelihood function; 2) the second part of the work is focused on the validation of the developed models by providing a detailed comparison with the present SoA, based on the confidence intervals of the parameters that maximize the quasi-likelihood function; 3) in the third-last part, the authors illustrate an application of the developed methodology to the disposable data coming from Tuscany Region1 to show its capability compared to the available SoA; 4) finally, discussion and conclusions are provided with improvements and future development. All the numerical analyses and computations described in this work have been implemented with the Matlab software package.

1.2. Structure and main ideas

The Generalized Linear Model (GLM) in the last few years has become a significant reality in claims reserve analysis [7,8,14,15,25]. The present work supplies the context of their forecast and applicability. The authors found the necessity, at first, to better describe and elucidate the mathematical background inside the subject and to assemble a coherent and well-defined model. In this context, the authors contributed to creating a new framework that could provide adequate and complete results in claims reserve.

The basis to estimate the claims reserve is through the construction of the Run-Off triangles [20] that estimate the total paid amount of claims incurred that are not yet reported to the insurance agency. A common methodology to estimate claims reserve is the Chain-Ladder (CL) [12,21,26]. The CL method aims to approximate the value of claims reserves that have been incurred but not yet reported and to predict the ultimate loss amounts in general insurance. The fundamental assumption on which the CL is based is that the patterns observed in historical loss development must hold for future loss development. The previous assumption represents one of the limitations of the CL and other deterministic methods because they do not extrapolate the tail factors [8] to increase the prediction performance for the paid amount. Several modifications of the CL have been proposed in the literature to improve the prediction of the claims reserve [11,23], and various stochastic claim reserving models have been developed for the general insurance to extrapolate the tail factors.

The employment of the stochastic models replicates the conventional reserve estimates produced by the CL methods, including parametric curves and smoothing techniques to reshape the runoff triangles, enabling the estimation of the tail factors through extrapolation. Stochastic models offer a significant benefit in providing precision measures for claims reserve, particularly emphasising the mean squared error of the prediction. However, understanding the stochastic distribution function of the claims reserve outcome is crucial to provide good prediction performances [19]. In this regard, the GLM has become, in the last years, a robust stochastic methodology for predicting claims reserve [12,17]. The GLM represents a wide range of models that allow the definition and estimation of parameters by maximising the likelihood function [17]. In this framework, the likelihood function is employed to evaluate the goodness of fit for a distribution that belongs to the exponential family of random variables. To define the likelihood function, one must specify the analytical form of the response variable distribution. In the works of [27,28], the last assumption has been overcome by introducing the quasi-likelihood function, which only needs the definition of the mean-variance relation of the response variable without any assumption on the shape of its distribution. Following the work of [22], the authors assume a logarithmic mean-variance relation, which provides the well-known over-dispersed Poisson (ODP) distribution model. The maximisation of the ODP quasi-likelihood function is performed using the Gauss-Newton (GN) algorithm [9,22]. This methodology has several benefits, but at the same time, it presents some crucial weaknesses that need to be considered and addressed.

Concerning the SoA, we propose a methodology to compute parameters in the framework of the GLM, which allows making good predictions for expectation value and errors without increased computational costs. Compared with other methods proposed in the SoA, our approach allows forecasting the claims reserve values and the corresponding standard errors by optimising two different residues. Consequently, it gives a better overview of the claims reserve in the healthcare sector, where the strong aleatory and peculiarity of claims’ lifespan need particular attention. The authors performed specific tests on different data-sets to validate the developed model. Firstly, the authors conducted a validation of the developed model to estimate the claims reserve to the data-sets present in the SoA, particularly from [22] and [23]. These validation tests help to understand if the developed model outperforms the ones in the literature. Subsequently, the developed model has been tested in the health insurance claims database of the Tuscany Region (Italy), which covers the period from 2010 to 2021. Dealing with this kind of database in the healthcare system is extremely rare and represents a flaw in estimating the claims reserve.

2. Theoretical background

In this Section, the authors present the theoretical background to estimate the claims reserve. The first Section (Section 2.1) describes the CL approach and how to construct the Run-Off triangles that represent the applied basis of the developed methodology. The second Subsection (Section 2.2) provides an introduction to GLM to understand the context in which our methodology is applied. The third Subsection (Section 2.3) introduces the definition of the quasi-likelihood function to make predictions of the response variable defined in the GLM framework. The fourth Subsection (Section 2.4) describes the implementation of the Gauss-Newton (GN) algorithm that maximises the quasi-likelihood function and its limitation. The fifth Subsection (Section 2.5) depicts the deployment of the genetic algorithm (GA), which overcomes the limitation problems of the GN algorithm, to provide consistent results and predictions. Finally, the sixth Subsection (Section 2.6) provides a workflow on how to estimate the errors of the quasi-likelihood function.

2.1. Run-Off triangles and Chain-Ladder approach

The characteristic feature of claims is their recognizable time lag from the point of incidence until they are notified to the insurance agency. As the reported claims are documented yearly, the insurance agency will face a significant amount of claims in a given calendar year, where each of them has a typical time delay. In this context, the Run-Off triangles [20] estimate the total paid amount of claims incurred during a period (e.g., fiscal year) that are not yet reported, and they help insurance agencies to estimate the amount of money that they will need to set aside to cover future claims that have already been incurred but have not been reported. For this reason, a reserve must be set aside for this purpose.

The Run-Off triangle is divided into cells where the row corresponds to the $i^{th}$ “incurred” year with $i\in [1,I]$, i.e. the year for which the claims happen, and the column corresponds to the $j^{th}$ “progression” year with $j\in [0,J]$, i.e. the delayed year of the payments from when the claims have been incurred. Each cell of the Run-Off triangle indicates the total paid amount $y_{i,j}$ for claims that have been incurred in the $i^{th}$ year and effectively paid in the $j^{th}$ year. In the following Table 1, the representation of the Run-Off triangle.

The most common approach used to calculate the claims reserve is the CL approach. Considering the definition of the Run-Off triangle represented in Table 1, and according to [22] and [20], the first step of the CL approach is the definition of the Run-Off triangle of the cumulative payments, which is shown in the following equation 1:

$$C_{i,j}=\sum _{k=0}^j{y}_{i,k}$$

(1)

where the variable $C_{i,j}$ represents the cumulative payment concerning the incurred-year i relative to the delayed-year j. The idea behind the CL method is that the cumulative payments of two adjacent anti-duration years, i.e. the $y_{i,j-1}$ and $y_{i,j}$ year, are proportional to each other, barring an erratic component $ϵ$ with zero expectation value. Another vital assumption of the CL is that the system is not susceptible to structural changes over time. Therefore, the estimate of the ultimate cumulative paid ${\widehat{C}}_{i,J}$ is obtained through the dilation factors ${\widehat{f}}_j^{\left(t\right)}$ (called link-ratio), according to equation 2:

$${\widehat{C}}_{i,J}^{\left(t\right)}=C_{i,t-i}\prod _{j=t-i}^{J-1}{\widehat{f}}_j^{\left(t\right)}$$

(2)

where the index $t\in [1,I]$ indicates the generic budget year, while the term ${\widehat{f}}_j^{\left(t\right)}$ represents the link-ratio which are defined through the following equation 3:

$${\widehat{f}}_j^{\left(t\right)}={\displaystyle \frac{{\sum }_{k=1}^{t-j-1}{C}_{k,j+1}}{{\sum }_{k=1}^{t-j-1}{C}_{k,j}}}$$

(3)

Hence, given the equations 2 and 3, the CL method is able to estimate the claims reserve ${\widehat{R}}_i^{\left(t\right)}$ given by equation 4:

$${\widehat{R}}_i^{\left(t\right)}={\widehat{C}}_{i,J}^{\left(t\right)}-C_{i,t-1}$$

(4)

In addition, in order to estimate the claims reserve per year, the stochastic approach allows the calculation of the complete claims reserve, which predicts the total payment amount for all the incurred and delayed years. In the following Sections, the authors will show the procedure to develop the stochastic model with the related theoretical approach to estimate the entire claims reserve, including the associated errors.

2.2. Generalized Linear Model

The GLM was introduced in the following work of [17] as a methodology to unify and generalize various statistical methods such as linear regression, logistic regression, etc. In ordinary linear regression, given a random variable Y, called response variable, one can assume that it is defined as a linear combination of m observations $y_i$, called explicative variables or predictors that are defined by equation 5:

$$Y=p(y_1,\dots ,y_m)=\sum _{i=1}^m{\gamma }_i{y}_i+\epsilon$$

(5)

where $\epsilon$ is a small variations with expectations value $E\left[\epsilon \right]=0$ and variance var$\left[\epsilon \right]={\sigma }^2$, while ${\gamma }_i$ are the coefficients of the linear combinations of predictors. In this framework, the estimation of Y became an algebraic problem of inverting matrices, but despite the simplicity of this methodology, this method has strong limitations, specifically when:

• the function p is not linear;
• the explicative variables $y_i$ don’t take values from the interval $(-\infty ,\infty )$;
• the variance ${\sigma }^2$ is not constant;
• it is not possible to assume the response variable Y following the normal distribution.

The GLM allows to overcome previous limitations by expanding the sets of applicability of the standard Linear Model and became a common tool in statistical analysis and prediction theory. Indeed, the GLM need a unique starting hypothesis about the behaviour of the probability distribution function (PDF) of Y, which is only required to belong to the class of exponential PDF, which includes: Normal, Poisson, Binomial, Gamma and Inverse Gaussian distribution. Moreover, if more response variables (Y) exist, they must be stochastically independent. Therefore, in general terms, the authors consider a response variable that follows a PDF parameterized by $\beta$ as defined in equation 6:

$$f(y;\beta ,\varphi )=e^{\alpha \left(\varphi \right)\left[y\beta -g\left(\beta \right)+h\left(y\right)\right]+\alpha (\varphi ,y)}$$

(6)

where y stands for a single i-realization of Y. For the sake of simplicity, authors have introduced a one-index notation, different from the one used to represent the Run-Off triangles in Table 1. In this case, the authors moved from a matrix-like notation $y_{i,j}$ to a vector-like notation $y_i$. The variable $\varphi$ represents the nuisance parameter2 and $\alpha \left(\varphi \right)$ is a positive definite function. In order to define the momenta of the distributions, we introduce the likelihood $L\left(\beta \right|y)$ and the log-likelihood $\mathcal{L}\left(\beta \right|y)$ functions as shown in equation 7:

$$L\left(\beta \right|y)=ln\mathcal{L}\left(\beta \right|y)=ln\prod _{i=1}^{m}f(y_i;\beta ,\varphi )=\sum _{i=1}^{m}f(y_i;\beta ,\varphi )$$

(7)

Based on the results of [14], it is possible to derive two crucial properties of the log-likelihood function shown in equation 8:

$$\begin{array}{cc}\hfill E\left[{\displaystyle \frac{\partial \mathcal{L}}{\partial \beta }}\right]& =0\hfill \\ \hfill E\left[{\displaystyle \frac{{\partial }^2\mathcal{L}}{\partial {\beta }^2}}\right]& =-E{\left[{\displaystyle \frac{\partial \mathcal{L}}{\partial \beta }}\right]}^2\hfill \end{array}$$

(8)

which provides the statement shown in equation 9 for the generic exponential function (equation 7):

$$\begin{array}{cc}\hfill E\left[{\displaystyle \frac{\partial \mathcal{L}}{\partial \beta }}\right]& =E\left[\alpha \left(\varphi \right)[y-g^{\prime }\left(\beta \right)]\right]=0\hfill \\ \hfill & \Rightarrow E\left[y\right]\equiv \mu =g^{\prime }\left(\beta \right)\hfill \end{array}$$

(9)

together with the quadratic expectation value:

$$\begin{array}{cc}\hfill E\left[{\displaystyle \frac{{\partial }^2\mathcal{L}}{\partial {\beta }^2}}\right]& =E[-\alpha \left(\varphi \right)g^{″}\left(\beta \right)]\hfill \\ \hfill & \Rightarrow E\left[\alpha \left(\varphi \right)g^{″}\left(\beta \right)\right]=E\left[{\alpha }^2\left(\varphi \right){(y-\mu )}^2\right]\hfill \end{array}$$

(10)

and hence,

$$g^{″}\left(\beta \right)=\alpha \left(\varphi \right)E\left[{(y-\mu )}^2\right]$$

(11)

which finally provides the expression for the variance of y (equation 12):

$$\alpha \left(\varphi \right)\mathrm{var}\left(y\right)=g^{″}\left(\beta \right)=V\left(\mu \right)$$

(12)

It is easy to prove that for a response variable Y normally distributed, the GLM reduces to the ordinary linear model.

2.3. Quasi-likelihood function

To make predictions and estimations concerning the response variables Y, in the GLM frameworks, it is necessary at first to estimate parameters of the log-likelihood function, commonly with the least-squares method. This procedure has a substantial limitation: the necessity to make a forecast hypothesis on the specific shape of the distribution (equation 7). In the work of [27], a new methodology was proposed introducing the so-called quasi-likelihood function, which only needs to fix the relationship between the first two momenta of the distribution, so mean and variance. The mean-variance relation is entirely defined by the link function (g) as shown in equation 13.

$$\begin{array}{cc}\hfill \mathrm{var}\left(y\right)& ={\displaystyle \frac{1}{\alpha \left(\varphi \right)}}V\left(\mu \right)\hfill \\ \hfill & ={\displaystyle \frac{1}{\alpha \left(\varphi \right)}}{g}^{″}\left(\beta \right)\hfill \\ \hfill & ={\displaystyle \frac{1}{\alpha \left(\varphi \right)}}{g}^{″}\left({\left(g^{\prime }\right)}^{-1}\left(\mu \right)\right)\hfill \end{array}$$

(13)

In equation 13, the term ${\left(g^{\prime }\right)}^{-1}$ represents the inverse function of the first derivative of the link function (g) according to equation 9.

Following, we will assume, without any lack of generality, an exponential link function. Our proposed model robustly represents the well-known GLM with Over-Dispersed Poisson distribution (ODP). It is possible to prove that the overall function $\alpha \left(\varphi \right)$ is proportional to the inverse of the nuisance parameter multiplied by a proportional constant representing a weight parameter relative to the single realization of Y. For the sake of simplicity, we will assume all the weights equal to unity. Therefore, equation 13 can be written as indicated in equation 14:

$$\mathrm{var}\left(y\right)=\varphi g^{″}\left({\left(g^{\prime }\right)}^{-1}\left(\mu \right)\right)$$

(14)

Consequently, for an exponential link function, the mean-variance relation is given by equation 15:

$$E\left[y\right]=\mu ={\displaystyle \frac{1}{\varphi }}\mathrm{var}\left(y\right)=V\left(\mu \right)$$

(15)

Considering that $V\left(\mu \right)=\mu$ in equation 15, the nuisance parameter is estimated by the Pearson Correlation Coefficient (PCC):

$$\varphi \equiv \widehat{\varphi }={\displaystyle \frac{1}{m-n}}\sum _{i=1}^m{\displaystyle \frac{{(y_i-{\mu }_i)}^2}{V\left({\mu }_i\right)}}⟹\varphi ={\displaystyle \frac{1}{m-n}}\sum _{i=1}^m{\displaystyle \frac{{(y_i-{\mu }_i)}^2}{{\mu }_i}}$$

(16)

where m is the number of observable data, n is the number of parameters of the quasi-likelihood function, and the difference $m-n$, indeed, represents the degree of freedom (DOF) of the system.

At this point, given the i-realization of the response variable Y, the quasi-likelihood function $K(y_i,{\mu }_i)$ is defined in equation 17:

$${\displaystyle \frac{\partial K(y_i,{\mu }_i)}{\partial {\mu }_i}}={\displaystyle \frac{y_i-{\mu }_i}{\mathrm{var}\left(y_i\right)}}={\displaystyle \frac{y_i-{\mu }_i}{\varphi V\left({\mu }_i\right)}}$$

(17)

or equivalently employing its integral form (equation 18):

$$K(y_i,{\mu }_i)=\int {\displaystyle \frac{y_i-{\mu }_i^{\prime }}{\varphi V\left({\mu }_i^{\prime }\right)}}d{\mu }_i^{\prime }+c\left(y_i\right)$$

(18)

where the index i indicates the single observations. The variable $c\left(y_i\right)$ represents a function of $y_i$, which in the future will be omitted for simplicity. For a complete overview of the properties and theorems concerning the quasi-likelihood function, we refer to the work of [27]. Therefore, substituting equation 15 into equation 18 and integrating over the i-interval $[y_i,{\mu }_i]$, we obtain the following expression (equation 19) for the quasi-likelihood function of the ODP model:

$$\begin{array}{cc}\hfill K(y_i,{\mu }_i)={\int }_{y_i}^{{\mu }_i}{\displaystyle \frac{y_i-{\mu }_i^{\prime }}{\varphi V\left({\mu }_i^{\prime }\right)}}d{\mu }_i^{\prime }& ={\int }_{y_i}^{{\mu }_i}{\displaystyle \frac{y_i-{\mu }_i^{\prime }}{\varphi {\mu }_i^{\prime }}}d{\mu }_i^{\prime }\hfill \\ \hfill & ={\displaystyle \frac{1}{\varphi }}\left[y_{i}ln{\displaystyle \frac{{\mu }_i}{y_i}}-{\mu }_i+y_i\right]\hfill \end{array}$$

(19)

Finally, summing over the complete set of the Y realization, it is possible to obtain the complete quasi-likelihood function in terms of y, $\beta$, and $\varphi$:

$$K(y,\beta ,\varphi )={\displaystyle \frac{1}{\varphi }}\sum _{i=1}^m\left[y_{i}ln{\displaystyle \frac{{\mu }_i}{y_i}}-{\mu }_i+y_i\right]$$

(20)

At this point, in the quasi-likelihood analysis, to obtain correct predictions for the response variable, it is necessary to estimate the parameters which maximize the function $K(y_i,{\mu }_i)$. In this regard, the non-linear link function requires an iterative algorithm to solve the problem. The adopted solution to estimate the parameters which maximize the function $K(y_i,{\mu }_i)$ is the Gauss-Newton algorithm (GN). The following Section will present the GN iterative methodology to estimate the parameters that maximize the $K(y_i,{\mu }_i)$ function. For clarity, the variable ${\mu }_i$ will be indicated as ${\mu }_i\left(\beta \right)$ because it represents the model that contains the parameters $\beta$ that maximize the quasi-likelihood function. Additionally, the details of its convergence problems and a solution to this problem will be proposed for a better choice of the initial guesses. In the following Sections, the authors will employ the vectorial notation for the $\beta$ parameters for completeness.

2.4. Gauss-Newton algorithm

Taken $\mathit{\beta }=({\beta }_1,\cdots ,{\beta }_n)$ the vector of parameters from which depend the response variable, to determine the ${\mathit{\beta }}^{*}$ vector that maximizes the quasi-likelihood function (see equation 19) we need to solve the following system of non-linear equation 21:

$${\displaystyle \frac{\partial K(y_i,{\mu }_i\left(\mathit{\beta }\right))}{\partial {\beta }_k}}=0$$

(21)

An iterative algorithm is required to locate the $k^{th}$ parameters of the $\mathit{\beta }$ vector, where $k\in [1,n]$, that solve the non-linear equation system shown in equation 21. Therefore, the GN algorithm has been deployed based on the guidelines provided in the work of [27]. Given a set of generic function (called residues) ${\left\{r_i\left(\mathit{\beta }\right)\right\}}_{i=1}^m$, the GN locates the set of ${\beta }_k$ parameters that minimize the function $S\left(\mathit{\beta }\right)$ defined in equation 22:

$$S\left(\mathit{\beta }\right)=\sum _{i=1}^m{r}_i^2\left(\mathit{\beta }\right)$$

(22)

where the specific choice of the residues function $r_i\left(\mathit{\beta }\right)$ is crucial in the entire process of parameter estimation. Following, the authors tested two different $r_i\left(\mathit{\beta }\right)$ functions, which we will define later, providing different results. Starting from the quasi-likelihood function $K(y_i,{\mu }_i\left(\mathit{\beta }\right))$, the Taylor expansion around the parameters $\mathit{\beta }$ is needed to minimize the function $S\mathit{\left(}\mathit{\beta }\right)$. Considering an expansion step $\delta$ given by equation 23:

$$\delta ={\mathit{\beta }}^{*}-\mathit{\beta }$$

(23)

The first order Taylor expansion in $\delta$ can be written in equation 24:

$$K(y_i,{\mu }_i\left({\mathit{\beta }}^{*}\right))=K(y_i,{\mu }_i\left(\mathit{\beta }\right))+\frac{\partial K(y_i,{\mu }_i\left(\mathit{\beta }\right))}{\partial \mathit{\beta }}{|}_{\mathit{\beta }}\delta +o\left({\delta }^2\right)$$

(24)

It is important to note that, because the function $K(y_i,{\mu }_i\left(\mathit{\beta }\right))$ is maximized by ${\mathit{\beta }}^{*}$, its first derivative in equation 24 is null when $\mathit{\beta }={\mathit{\beta }}^{*}$, therefore a maximization procedure must be developed to locate the ${\delta }^{*}$ parameters. The GN works directly on the function $S\left(\mathit{\beta }\right)$ defined in equation 22, which has no triviality problem. Subsequently, it remains to define a consistent definition of the function $r_i\left(\mathit{\beta }\right)$ to start the optimization procedure properly and minimize the $S\left(\mathit{\beta }\right)$ function through the GN methodology. In this Section, we propose two different definitions of residues $r_i\left(\mathit{\beta }\right)$, which allows to work separately on forecasting the claims reserve and their corresponding standard errors. The two definitions are consistent with providing the best estimation for the reserves and trying to reduce the related CI.

The first formulation of the residues, as characterised by [22], concerns the direct optimization of the parameters inside the quasi-likelihood $K(y_i,{\mu }_i\left(\mathit{\beta }\right))$, and it is given by equation 25:

$$r_i\left(\mathit{\beta }\right)=\frac{y_i}{\varphi }-K(y_i,{\mu }_i\left(\mathit{\beta }\right))⟹r_i\left(\mathit{\beta }\right)={\displaystyle \frac{1}{\varphi }}\left[{\mu }_i\left(\mathit{\beta }\right)-y_{i}ln{\displaystyle \frac{{\mu }_i\left(\mathit{\beta }\right)}{y_i}}\right]$$

(25)

The second formulation of the residues, developed by the authors, concerns the optimization of the dispersion parameters ($\varphi$ - see equation 16), and it is given by the equation 26:

$$r_i\left(\mathit{\beta }\right)={\displaystyle \frac{y_i-{\mu }_i\left(\mathit{\beta }\right)}{\sqrt{{\mu }_i\left(\mathit{\beta }\right)}}}$$

(26)

Regarding the following steps, the authors will not use one of the previous definitions of residues but will keep the discussion in general terms. The comparison of the two specific formulations of $r_i\left(\mathit{\beta }\right)$ (equations 25, and 26) will be discussed in Section 3 and Section 4, where the results related to comparison with the SoA, and the data from the Tuscany Region are presented, respectively. The function $S\left(\mathit{\beta }\right)$ does not present triviality problems compared to the quasi-likelihood function. Therefore, the first order Taylor expansion of equation 22 can be performed in $\delta$ (see equation 23) around the initial point $\mathit{\beta }$ using a vectorial formalism provides the equation 27:

$$\begin{array}{cc}\hfill S\left({\mathit{\beta }}^{*}\right)& =S\left(\mathit{\beta }\right)+{\displaystyle \frac{\partial S\left(\mathit{\beta }\right)}{\partial \mathit{\beta }}}{|}_{\mathit{\beta }}\delta +o\left({\delta }^2\right)\hfill \\ \hfill & =S\left(\mathit{\beta }\right)+2{\mathit{r}}^T\left(\mathit{\beta }\right)·{\displaystyle \frac{\partial \mathit{r}\left(\mathit{\beta }\right)}{\partial \mathit{\beta }}}{|}_{\mathit{\beta }}\delta +o\left({\delta }^2\right)\hfill \end{array}$$

(27)

The first derivative in equation 27 represents the Jacobian matrix of the residues computed in the initial point $\mathit{\beta }$. More precisely, the first derivative appearing in equation 27 depends on the definition of residues (see equations 25 and 26). Therefore the single elements of the Jacobian matrix can be defined by the following equations (28-29):

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial r_i\left(\mathit{\beta }\right)}{\partial {\mathit{\beta }}_k}}{|}_{\mathit{\beta }}& =-{\displaystyle \frac{1}{\varphi }}{\displaystyle \frac{\partial {\mu }_i\left(\mathit{\beta }\right)}{\partial {\mathit{\beta }}_k}}{|}_{\mathit{\beta }}{\displaystyle \frac{y_i-{\mu }_i\left(\mathit{\beta }\right)}{{\mu }_i\left(\mathit{\beta }\right)}}\equiv J_{ik}\left(\mathit{\beta }\right)(\mathrm{see}\mathrm{equation}25)\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial r_i\left(\mathit{\beta }\right)}{\partial {\mathit{\beta }}_k}}{|}_{\mathit{\beta }}& =-{\displaystyle \frac{1}{2}}{\displaystyle \frac{\partial {\mu }_i\left(\mathit{\beta }\right)}{\partial {\mathit{\beta }}_k}}{|}_{\mathit{\beta }}{\displaystyle \frac{y_i+{\mu }_i\left(\mathit{\beta }\right)}{{\mu }_i{\left(\mathit{\beta }\right)}^{\frac{3}{2}}}}\equiv J_{ik}\left(\mathit{\beta }\right)(\mathrm{see}\mathrm{equation}26)\hfill \end{array}$$

(29)

The term appearing in the left-hand side of the equations 28 and 29 represents the derivative of the $i^{th}$-residue $r_i\left(\mathit{\beta }\right)$ space with respect to every ${\mathit{\beta }}_k$-parameters. The variable $J_{ik}\left(\mathit{\beta }\right)$ appearing in the right-hand side of the equations 28 and 29 represents the i-row and k-column of the Jacobian matrix defined in equation 27, thus the Jacobian matrix has a dimension of $m\times n$ where m is the number of residues (i.e. the observations - see equation 5) and n represents the number of parameters (i.e. the length of the $\mathit{\beta }$ vector). At this step, the previous Taylor expansion shown in equation 27 can be rewritten as follows (equation 30):

$$\begin{array}{cc}\hfill S\left({\mathit{\beta }}^{*}\right)& =S\left(\mathit{\beta }\right)+2{\mathit{r}}^T\left(\mathit{\beta }\right)·\mathit{J}\left(\mathit{\beta }\right)·\delta +o\left({\delta }^2\right)\hfill \\ \hfill & ={\mathit{r}}^T\left(\mathit{\beta }\right)·\mathit{r}\left(\mathit{\beta }\right)+2{\mathit{r}}^T\left(\mathit{\beta }\right)·\mathit{J}\left(\mathit{\beta }\right)·\delta +o\left({\delta }^2\right)\hfill \end{array}$$

(30)

The term $\mathit{J}\left(\mathit{\beta }\right)$ appearing in equation 30, represents the Jacobian matrix in a vectorial formalism computed in the initial point $\mathit{\beta }$. The variable $S\left(\mathit{\beta }\right)$ in equation 30 is expressed by performing the dot product between the residues $\mathit{r}\left(\mathit{\beta }\right)$ and its vectorial transpose. For the sake of simplicity, the variables $\mathit{J}\left(\mathit{\beta }\right)$ and $\mathit{r}\left(\mathit{\beta }\right)$ will be indicated as $\mathit{J}$ and $\mathit{r}$, respectively. At this point, it is possible to compute the gradient of the $S\left(\mathit{\beta }\right)$ function in the stationary point ${\mathit{\beta }}^{*}$ as appearing in equation 30:

$${\displaystyle \frac{\partial S\left(\mathit{\beta }\right)}{\partial \mathit{\beta }}}{|}_{{\mathit{\beta }}^{*}}=2{\mathit{J}}^T\mathit{r}+2{\mathit{J}}^T\mathit{J}\delta +2{\mathit{r}}^T{\displaystyle \frac{\partial \mathit{J}}{\partial \mathit{\beta }}}{|}_{\mathit{\beta }}\delta$$

(31)

Since the gradient of $S\left(\mathit{\beta }\right)$ is computed in the stationary point ${\mathit{\beta }}^{*}$, the derivatives of this gradient are null. Additionally, equation 31 can be further simplified by neglecting the term containing the derivative of the Jacobian matrix. Therefore, equation 31 can be rewritten as follows:

$${\displaystyle \frac{\partial S\left(\mathit{\beta }\right)}{\partial \mathit{\beta }}}{|}_{{\mathit{\beta }}^{*}}=2{\mathit{J}}^T\mathit{r}+2{\mathit{J}}^T\mathit{J}\delta =0$$

(32)

Therefore, the $\delta$ can be re-calculated as function of ${\mathit{\beta }}^{*}$ (see equation 23) as follow:

$$\delta =-{\left({\mathit{J}}^T\mathit{J}\right)}^{-1}{\mathit{J}}^T\mathit{r}⟹{\mathit{\beta }}^{*}=\mathit{\beta }-{\left({\mathit{J}}^T\mathit{J}\right)}^{-1}{\mathit{J}}^T\mathit{r}$$

(33)

In practice, the GN algorithm computes a ${\mathit{\beta }}^{\left(s\right)}$ value at every iterative step (s) starting from an initial guess of ${\mathit{\beta }}^{\left(0\right)}$. Therefore, the ${\delta }^{\left(s\right)}$ value is calculated using equation 34:

$${\delta }^{\left(s\right)}=-{\left({\mathit{J}}^T\left({\mathit{\beta }}^{\left(s\right)}\right)·\mathit{J}\left({\mathit{\beta }}^{\left(s\right)}\right)\right)}^{-1}{\mathit{J}}^T\left({\mathit{\beta }}^{\left(s\right)}\right)·\mathit{r}\left({\mathit{\beta }}^{\left(s\right)}\right)$$

(34)

The iterative procedure described in equation 34 is carried on until the following condition is verified (equation 35):

$$∥{\delta }^{\left(s\right)}∥\equiv ∥{\mathit{\beta }}^{(s+1)}-{\mathit{\beta }}^{\left(s\right)}∥<ϵ$$

(35)

where the variable $ϵ$ represents an arbitrary user-defined tolerance for which the GN algorithm must reach to stop the iterations. At this point, the GN algorithm has found the $\mathit{\beta }$ parameters, which minimize the function $S\left(\mathit{\beta }\right)$ and therefore, the final values of $\mathit{\beta }$ is given as follows (equation 36):

$${\mathit{\beta }}^{*}\equiv {\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{\left(s\right)}-{\left({\mathit{J}}^T\left({\mathit{\beta }}^{\left(s\right)}\right)·\mathit{J}\left({\mathit{\beta }}^{\left(s\right)}\right)\right)}^{-1}{\mathit{J}}^T\left({\mathit{\beta }}^{\left(s\right)}\right)·\mathit{r}\left({\mathit{\beta }}^{\left(s\right)}\right)$$

(36)

The described procedure represents an iterative methodology to locate the stationary point that maximizes the quasi-likelihood function $K(y_i,{\mu }_i\left(\mathit{\beta }\right))$, and it can be iterated until the $∥\delta ∥$ reaches an arbitrary tolerance value ($ϵ$ = 10⁻¹⁵) defined by the user. For each step of the iterative method, the $\mathit{\beta }$ is updated closer to the value that minimizes the $S\left(\mathit{\beta }\right)$ function. The choice of a good initial guess for $\mathit{\beta }$ is a crucial point for the iterative GN method to have success in locating the global minimum of the $S\left(\mathit{\beta }\right)$ function. However, estimating a good initial guess of the multi-parametric function $S\left(\mathit{\beta }\right)$ to force the convergence of the GN method to a global minimum requires additional concerns and tools that must be addressed.

2.4.1. The Hessian modification

The convergence of the GN algorithm (equation 36) is not guaranteed for the residues defined in equations 25 and 26. The main reason that the GN algorithm fails to reach the convergence is due to the approximation introduced in equation 32 where the derivative of the Jacobian matrix ($\mathit{J}$) with respect to $\mathit{\beta }$ is ignored. In order to be able to apply the equation 36, the following approximation must hold:

$$\left|{\mathit{J}}^T\mathit{J}\right|\gg \left|{\mathit{r}}^T{\displaystyle \frac{\partial \mathit{J}}{\partial \mathit{\beta }}}\right|$$

(37)

In order to expect the convergence of the GN algorithm using the equation 36, the approximation highlighted in equation 37 must be valid in one of the two cases: 1) the residues $\mathit{r}\left(\mathit{\beta }\right)$ are small compared to the derivatives of the Jacobian ($\partial \mathit{J}/\partial \mathit{\beta }$) around the stationary point ${\mathit{\beta }}^{*}$; 2) the residues $\mathit{r}\left(\mathit{\beta }\right)$ are blandly non-linear in such a way that the derivatives of the Jacobian are negligible [13]. Therefore, in the case for which the approximation 37 cannot hold, the key equation of the GN algorithm (equation 36) must be modified. Starting from equation 31, it can be rearranged in the following equation:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial S\left(\mathit{\beta }\right)}{\partial \mathit{\beta }}}{|}_{{\mathit{\beta }}^{*}}& =2{\mathit{J}}^T\mathit{r}+2\left[{\mathit{J}}^T\mathit{J}+2{\mathit{r}}^T{\displaystyle \frac{\partial \mathit{J}}{\partial \mathit{\beta }}}{|}_{\mathit{\beta }}\right]·\delta \hfill \\ \hfill \nabla S\left({\mathit{\beta }}^{*}\right)& =\nabla S\left(\mathit{\beta }\right)+\mathit{H}\left(S\left(\mathit{\beta }\right)\right)·\delta \hfill \end{array}$$

(38)

The variable $\nabla S\left(\mathit{\beta }\right)$ in equation 38 represents the gradient of the function $S\left(\mathit{\beta }\right)$ (see equation 22) computed in the initial point $\mathit{\beta }$, the term $\mathit{H}\left(S\left(\mathit{\beta }\right)\right)$ indicates the Hessian of the function $S\left(\mathit{\beta }\right)$ computed in the initial point $\mathit{\beta }$. The Hessian is a matrix with a dimension of $n\times n$ where n n represents the number of parameters (i.e. the length of the $\mathit{\beta }$ vector), where each element of this matrix contains the second derivatives of the function $S\left(\mathit{\beta }\right)$. For simplicity the gradient $\nabla S\left(\mathit{\beta }\right)$ will be indicated as $\mathit{g}$, and the Hessian $\mathit{H}\left(S\left(\mathit{\beta }\right)\right)$ will be indicated as $\mathit{H}$. Since the gradient $\nabla S\left({\mathit{\beta }}^{*}\right)$ is computed in the stationary point ${\mathit{\beta }}^{*}$, the derivatives of this gradient are null. Therefore, equation 38 can be rewritten as follows:

$$\begin{array}{cc}\hfill \nabla S\left({\mathit{\beta }}^{*}\right)& =\mathit{g}+\mathit{H}·\delta =0⟹\hfill \\ \hfill \delta & =-{\mathit{H}}^{-1}·\mathit{g}\hfill \end{array}$$

(39)

Therefore, the $\delta$ can be re-calculated as function of ${\mathit{\beta }}^{*}$ (see equation 23) as follow:

$${\mathit{\beta }}^{*}=\mathit{\beta }-{\mathit{H}}^{-1}·\mathit{g}$$

(40)

In this case, the GN algorithm computes a ${\mathit{\beta }}^{\left(s\right)}$ value at every iterative step (s) starting from an initial guess of ${\mathit{\beta }}^{\left(0\right)}$. Therefore, the ${\delta }^{\left(s\right)}$ value is calculated using the following equation 41:

$${\delta }^{\left(s\right)}=-{\left(\mathit{H}\left({\mathit{\beta }}^{\left(s\right)}\right)\right)}^{-1}·\mathit{g}\left({\mathit{\beta }}^{\left(s\right)}\right)$$

(41)

The iterative procedure described in equation 41 is carried on until the condition $∥\delta ∥<ϵ$ is verified. The variable $ϵ$ represents an arbitrary user-defined tolerance for which the GN algorithm must reach to stop the iterations. At this point, the GN algorithm has found the $\mathit{\beta }$ parameters, which minimize the function $S\left(\mathit{\beta }\right)$ and therefore, the final values of $\mathit{\beta }$ is given as follows 36:

$${\mathit{\beta }}^{*}\equiv {\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{\left(s\right)}-{\left(\mathit{H}\left({\mathit{\beta }}^{\left(s\right)}\right)\right)}^{-1}·\mathit{g}\left({\mathit{\beta }}^{\left(s\right)}\right)$$

(42)

The described procedure represents an alternative iterative methodology to locate the stationary point (${\mathit{\beta }}^{*}$) that maximizes the quasi-likelihood function $K(y_i,{\mu }_i\left(\mathit{\beta }\right))$, and it can be adopted when the previous methodology (see Section 2.4) fails completely to locate the stationary point.

2.4.2. Limits of the GN algorithm

The convergence of the GN algorithm using the Jacobian matrix (see Section 2.4) and its Hessian modification (see Section 2.4.1) is not always guaranteed, as proved and discussed by [13] and [16]. Even local convergence could fail for the following reasons: 1) the initial guess ${\mathit{\beta }}^{\left(0\right)}$ is far from the stationary point ${\mathit{\beta }}^{*}$; 2) the ${\mathit{J}}^T\mathit{J}$ matrix described in equation 33 is ill-conditioned; the Hessian matrix $\mathit{H}$ in equation 40 is ill-conditioned. Therefore, in order to perform claims reserve in this framework, it is necessary to address and cure the problem of GN convergence and provide a method to ensure efficiency in finding the global minimum by the algorithm. Local convergence allows predicting acceptable results for the claims reserve but could provide huge errors associated with the estimated $\mathit{\beta }$ parameters. Thus, providing an instrument that guarantees the GN algorithm’s global convergence is necessary to offer good results for claims reserve and associated errors without increasing computational costs. In this regard, the author aims to derive a good initial guess of the $\mathit{\beta }$ parameters to prevent divergence or local convergence of the GN algorithm. The genetic algorithm has been employed to help locate the good initial guess of the $\mathit{\beta }$ parameters to solve the convergence problem.

2.5. The genetic algorithm

The GA is a methodology for solving constrained and unconstrained optimization problems based on a natural selection process that mimics biological evolution. The GA alters a present population of the $\mathit{\beta }$ parameters by randomly picking and processing them as parents to create the next generation’s offspring. Over generations, the population evolves to the solution represented by a global minimum. The GA behaves differently from classical, derivative-based optimization algorithms, i.e. it can address issues that conventional optimization algorithms may not be able to process, such as discontinuous, non-differentiable, or highly non-linear objective functions. A detailed description and explanation of how the GA works is presented in the following references [4,5,10]. The algorithm, at each iteration, creates a sequence of new populations of individuals starting from the current generation.

The creation process is based on the following steps. Each population individual is associated with a score called "fitness" value, allowing the algorithm to select specific members. The fitness value is calculated by a fitness function represented by the residues defined in the equation 25 or equation 26. Individuals with high fitness are chosen as the elite population. Hence, they are selected as the elite population to create the next generation of individuals from which the children are produced by manipulating the $\beta$ parameters by combining them as an entry of a pair of parents-crossover or a single parent-mutation. In contrast, the algorithm excludes individuals with low fitness who are not participating in the selection process. Finally, the algorithm is stopped under arbitrary criteria that might result in reaching a certain deviation threshold in fitness values between the individuals, the maximum imposed number of generations, and the maximum number of stall generations corresponding to the number of generations where the best fitness value has not changed.

2.6. Expectation values and error estimation

Following the approach described in Section 2.4, the parameterization of the $\mathit{\beta }$ vector for the quasi-likelihood function defined in equation 18 must be carried out. The parameterization described by [22] has been followed to define the $\mathit{\beta }$ vector as shown in equation 43:

$$\mathit{\beta }=(c,a_1,\dots ,a_n,b_1,\dots ,b_m)$$

(43)

where the authors have used a different definition of the indexes n and m with respect to the Section 2.4, which is the one that will be used following. Therefore, the index $k\in [1,n+m+1]$ and it labels the k-parameters of the $\mathit{\beta }$ vector. The variable c represents the intercept, the variable $a_{n_i}$ (with $n_i\in [1,n]$) is associated with the $n_i^{th}$ rows of the Run-Off triangles, i.e. the incurred year, and the variable $b_{m_i}$ (with $m_i\in [1,m]$) corresponds to the $m_i^{th}$ columns of the Run-Off triangles, i.e. the progression year. In this case, the notation has the advantage that the subscript index (i) indicates the $i^{th}$ residue, which embeds the index couple $(n_i,m_i)$ where $n_i$ can be ≠ from $m_i$. In the ODP framework, the exponential function with the parameterization $\mathit{\beta }$, is represented by equation 44 and indicates the expectation value of the $y_i$ observable:

$$E\left[y_i\right]={\mu }_i\left(\mathit{\beta }\right)=\left\{\begin{array}{cc}{e}^c\hfill & \mathrm{if}{n}_i=m_i=1\hfill \\ e^{c+a_{n_i}}\hfill & \mathrm{if}{m}_i=1\mathrm{and}{n}_i\ne 1\hfill \\ e^{c+b_{m_i}}\hfill & \mathrm{if}{n}_i=1\mathrm{and}{m}_i\ne 1\hfill \\ e^{c+a_{n_i}+b_{m_i}}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(44)

Once the ODP framework has been chosen, the methodology discussed in Section 2.4 can be applied. The main goal is to find the best estimate of the $\mathit{\beta }$ parameters associated with the generic ${\mu }_i$ model (equation 44) to fill the Run-Off triangle and estimate the claims reserve jointly with their confidence interval. Therefore, it is necessary to find the parameters that maximise the quasi-likelihood function 18 by substituting it in the equation 44. In this manner, it is possible to obtain equation 45:

$$K(y_i,{\mu }_i\left(\mathit{\beta }\right))={\displaystyle \frac{1}{\varphi }}\left[y_i(c+a_{n_i}+b_{m_i})-y_{i}log{y}_i-e^{c+a_{n_i}+b_{m_i}}+y_i\right]$$

(45)

At this step, the function $S\left(\mathit{\beta }\right)$ defined in equation 22 has been built with the two residues $r_i\left(\mathit{\beta }\right)$ (equations 25-26) defined in the Section 2.4 and it depends directly on the parameters $\mathit{\beta }$. In the ODP framework, the residues are defined as the following equations (46-47):

$$\begin{array}{cc}\hfill r_i\left(\mathit{\beta }\right)={\displaystyle \frac{1}{\varphi }}\left[e^{c+a_{n_i}+b_{m_i}}+y_{i}log{y}_i-y_i(c+a_{n_i}+b_{m_i})\right]& (\mathrm{equation}25)\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill r_i\left(\mathit{\beta }\right)={\displaystyle \frac{y_i-e^{c+a_{n_i}+b_{m_i}}}{e^{\frac{1}{2}(c+a_{n_i}+b_{m_i})}}}& (\mathrm{equation}26)\hfill \end{array}$$

(47)

At this step, considering the total number of residues (N) which corresponds to the number of observations, the Jacobian matrix ($\mathit{J}$) to apply the GN algorithm (see equation 36) is defined as follows (equation 48):

$$\mathit{J}=\left(\begin{array}{ccccccc}{\displaystyle \frac{\partial r_1}{\partial c}}& {\displaystyle \frac{\partial r_1}{\partial a_1}}& \dots & {\displaystyle \frac{\partial r_1}{\partial a_n}}& {\displaystyle \frac{\partial r_1}{\partial b_1}}& \dots & {\displaystyle \frac{\partial r_1}{\partial b_m}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\displaystyle \frac{\partial r_i}{\partial c}}& {\displaystyle \frac{\partial r_i}{\partial a_1}}& \dots & {\displaystyle \frac{\partial r_i}{\partial a_n}}& {\displaystyle \frac{\partial r_i}{\partial b_1}}& \dots & {\displaystyle \frac{\partial r_i}{\partial b_m}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\displaystyle \frac{\partial r_N}{\partial c}}& {\displaystyle \frac{\partial r_N}{\partial a_1}}& \dots & {\displaystyle \frac{\partial r_N}{\partial a_n}}& {\displaystyle \frac{\partial r_N}{\partial b_1}}& \dots & {\displaystyle \frac{\partial r_N}{\partial b_m}}\end{array}\right)$$

(48)

Starting from the initial guess ${\mathit{\beta }}^{\left(0\right)}$ provided by the GA (Section 2.5), the GN algorithm updates the Jacobian matrix (equation 48) using one of the two residues defined in equations 46 or 47 to compute iteratively the parameters ${\mathit{\beta }}^{*}$ which minimize the function $S\left(\mathit{\beta }\right)$ through equation 36. During the iterative process using the Jacobian matrix, the relation 37 must hold to reach the convergence of the calculus. If equation 37 cannot hold the Hessian modification of the GN algorithm (Section 2.4.1) must be applied. In this case, the gradient ($\mathit{g}$) and the Hessian ($\mathit{H}$) of the function $S\left(\mathit{\beta }\right)$ are defined as follows (for simplicity, the variable $S\left(\mathit{\beta }\right)$ is denoted as S):

$$\mathit{g}={\left[\left({\displaystyle \frac{\partial S}{\partial c}}\right),\left({\displaystyle \frac{\partial S}{\partial a_1}}\right),\dots ,\left({\displaystyle \frac{\partial S}{\partial a_n}}\right),\left({\displaystyle \frac{\partial S}{\partial b_1}}\right),\dots ,\left({\displaystyle \frac{\partial S}{\partial b_m}}\right)\right]}^T$$

(49)

$$\mathit{H}=\left(\begin{array}{ccccccc}{\displaystyle \frac{{\partial }^{2}S}{\partial c^2}}& {\displaystyle \frac{{\partial }^{2}S}{\partial c\partial a_1}}& \dots & {\displaystyle \frac{{\partial }^{2}S}{\partial c\partial a_n}}& {\displaystyle \frac{{\partial }^{2}S}{\partial c\partial b_1}}& \dots & {\displaystyle \frac{{\partial }^{2}S}{\partial c\partial b_m}}\\ {\displaystyle \frac{{\partial }^{2}S}{\partial a_1\partial c}}& {\displaystyle \frac{{\partial }^{2}S}{\partial a_1^2}}& \dots & {\displaystyle \frac{{\partial }^{2}S}{\partial a_1\partial a_n}}& {\displaystyle \frac{{\partial }^{2}S}{\partial a_1\partial b_1}}& \dots & {\displaystyle \frac{{\partial }^{2}S}{\partial a_1\partial b_m}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\displaystyle \frac{{\partial }^{2}S}{\partial a_n\partial c}}& {\displaystyle \frac{{\partial }^{2}S}{\partial a_n\partial a_1}}& \dots & {\displaystyle \frac{{\partial }^{2}S}{\partial a_n^2}}& {\displaystyle \frac{{\partial }^{2}S}{\partial a_n\partial b_1}}& \dots & {\displaystyle \frac{{\partial }^{2}S}{\partial a_n\partial b_m}}\\ {\displaystyle \frac{{\partial }^{2}S}{\partial b_1\partial c}}& {\displaystyle \frac{{\partial }^{2}S}{\partial b_1\partial a_1}}& \dots & {\displaystyle \frac{{\partial }^{2}S}{\partial b_1\partial a_n}}& {\displaystyle \frac{{\partial }^{2}S}{\partial b_1^2}}& \dots & {\displaystyle \frac{{\partial }^{2}S}{\partial b_1\partial b_m}}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\displaystyle \frac{{\partial }^{2}S}{\partial b_m\partial c}}& {\displaystyle \frac{{\partial }^{2}S}{\partial b_m\partial a_1}}& \dots & {\displaystyle \frac{{\partial }^{2}S}{\partial b_m\partial a_n}}& {\displaystyle \frac{{\partial }^{2}S}{\partial b_m\partial b_1}}& \dots & {\displaystyle \frac{{\partial }^{2}S}{\partial b_m^2}}\end{array}\right)$$

(50)

At that point is possible to fill the entire Run-Off table with the complete set of estimations, given by equation 51:

$${\mu }_i\left({\mathit{\beta }}^{*}\right)=\left\{\begin{array}{cc}{e}^{c^{*}}\hfill & \mathrm{if}{n}_i=m_i=1\hfill \\ e^{c^{*}+a_{n_i}^{*}}\hfill & \mathrm{if}{m}_i=1\mathrm{and}{n}_i\ne 1\hfill \\ e^{c^{*}+b_{m_i}^{*}}\hfill & \mathrm{if}{n}_i=1\mathrm{and}{m}_i\ne 1\hfill \\ e^{c^{*}+a_{n_i}^{*}+b_{m_i}^{*}}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(51)

Concurrently with the estimation of the expectation values, the corresponding errors on the ${\mathit{\beta }}^{*}$ parameters have been computed starting from the Hessian matrix $\mathit{H}\left(K(y,{\mathit{\beta }}^{*},\varphi )\right)$ of the complete quasi-likelihood function defined in equation 20, where each element of the Hessian matrix ($H_{jk}$) is computed using equation 52:

$$H_{jk}={\displaystyle \frac{{\partial }^{2}K(y,{\mathit{\beta }}^{*},\varphi )}{\partial {\mathit{\beta }}_j\partial {\mathit{\beta }}_k}}$$

(52)

Therefore, the variance ($\delta {\mathit{\beta }}^2$) associated with each ${\mathit{\beta }}^{*}$ parameter is computed by performing the inverse of the Hessian matrix ($\mathit{H}$), and then taking the diagonal elements of the matrix changing the sign as shown in equation 53. For the sake of simplicity, the ${\mathit{\beta }}^{*}$ term is indicated as $\mathit{\beta }$, and the variable $\mathit{H}$ indicates the Hessian of the complete quasi-likelihood function $\mathit{H}\left(K(y,{\mathit{\beta }}^{*},\varphi )\right)$.

$$\delta {\mathit{\beta }}^2=\left(\delta c^2,\delta a_1^2,\dots ,\delta a_n^2,\delta b_1^2,\dots ,\delta b_m^2\right)=-\mathrm{diag}\left({\mathit{H}}^{-1}\right)$$

(53)

In order to provide the complete estimation of the standard errors for the expectation values ${\mu }_i\left(\mathit{\beta }\right)$, the authors applied the theory of error propagation. Therefore, they computed the differential of the generic expectation values ${\mu }_i\left(\mathit{\beta }\right)$ (see equation 51), which is given by equation 54:

$$d{\mu }_i=d\left(e^{c+a_{n_i}+b_{m_i}}\right)=e^{c+a_{n_i}+b_{m_i}}(dc+d{a}_{n_i}+d{b}_{m_i})$$

(54)

Hence, substituting the parameters errors $\delta c$, $\delta a_{n_i}$ and $\delta b_{m_i}$ with the corresponding differentials in equation 54:

$$\delta {\mu }_i=e^{c+a_{n_i}+b_{m_i}}(\delta c+\delta a_{n_i}+\delta b_{m_i})$$

(55)

The variance $\delta {\mu }_i$ of the expectation values is given by the squared sum of the single term of the previous equation 55:

$$\delta {\mu }_i^2=e^{2(c+a_{n_i}+b_{m_i})}(\delta c^2+\delta a_{n_i}^2+\delta b_{m_i}^2+2\delta c\delta a_{n_i}+2\delta c\delta b_{m_i}+2\delta a_{n_i}\delta b_{m_i})$$

(56)

Where the mixed terms are negligible for independent variables. In this framework the standard errors of the expectation value $\Delta {\mu }_i$ is defined as follows 57:

$$\Delta {\mu }_i=e^{c+a_{n_i}+b_{m_i}}\sqrt{d{c}^2+d{a}_{n_i}^2+d{b}_{m_i}^2}$$

(57)

This Section completes the treatment of the GLM methods using the GN algorithm to estimate parameters with an initial guess provided via the GA algorithm. Starting from the SoA approach, the procedure developed in these Sections (Section 2.3, Section 2.4 and Section 2.5) is used by the authors to provide results in the context of the healthcare sector, and specifically for the Tuscany Region case study. The aim is to provide results in the form given by the theoretical background of Section 2.2, composed by expectation value and corresponding standard error. The authors are now able to present, in the next Section, the validation of the developed approach by comparison with the available SoA.

3. Comparison with SoA Models

To validate the proposed methodology, the authors performed a comparison with two works present in the SoA. The first comparison concerns the work of [22], which primarily inspired the authors because it presents a similar methodology, presenting a sweetened data-set from a different context from the healthcare (i.e., RC cars in Italy). The data-set of [22] provides a complete calculation of the expectation values with the corresponding errors to the estimated parameters ($\mathit{\beta }$). The second comparison concerns the work of [23], which is a good data-set that contains predictions coming from a robustification of the CL method with highly homogeneous data. The aim is to apply the methodology developed in the previous Sections (i.e. Section 2.4, Section 2.5, and Section 2.6) and improve the results compared with the ones in the SoA. Following, the authors show two different figures for each of the mentioned works to compare the results coming from the different choices of the residues, defined respectively in equations 46 and . The key performance indicators (KPI) used to compare our methodology with the ones proposed in the literature are the function $S\left(\mathit{\beta }\right)$ computed as shown in equation 22 using the values of the $\mathit{\beta }$ vector that minimize the residues after applying the GN algorithm (${\mathit{\beta }}^{*}$), the dispersion parameter ($\varphi$) computed as shown in equation 16, and the claims reserve amount ($R_{tot}$) that must be set aside including the corresponding errors of the expectation value.

The developed stochastic model allows to compute claims reserve differently with respect to the CL method described in Section 2.1. After computing the ${\mathit{\beta }}^{*}$ parameters of the stationary point through the GN algorithm, the expectation values shown in equation 51 can be calculated and they correspond to every single cell of the Run-Off triangle. Considering the two-indices notation adopted in the Run-Off triangle of Table 1, the developed stochastic approach is able to calculate the complete set of expectation values in order to estimate the entire claims reserve. For clarity, the expectation values are indicated in the following Table 2.

The coloured cells in "light blue" as shown in Table 2 represent the expectation values ${\mu }_{i,j}$ with $i+j\le I$ coming from the application of the GN algorithm. In other words, they correspond to the generic model ${\mu }_i\left({\mathit{\beta }}^{*}\right)$ indicated in Section 2.6, for which they have been used to construct the residues $r_i\left(\mathit{\beta }\right)$ in the equation 46 and 47. Additionally, these expectation values (${\mu }_{i,j}$) are important to understand the dispersion of the observable data ($y_{i,j}$ in Table 1) and therefore compute the dispersion parameter ($\varphi$). In contrast, the coloured cells in "light orange" as shown in Table 2 represent the expectation values ${\mu }_{i,j}$ with $i+j>I$ that are estimated using equation 51, and they are being used to predict the claims reserve for the generic budget year ($t>I$) through equation 58:

$$R^{\left(t\right)}=\sum _{i+j=t}{\mu }_{i,j}\mathrm{with}t\in [I+1,I+J]$$

(58)

In order to predict the entire claims reserve $R_{tot}$, i.e. the total amount to be set aside, is necessary to sum over every budget year $t\in [I+1,I+J]$, through equation 59:

$$R_{tot}=\sum _{t=I+1}^{I+J}{R}^{\left(t\right)}$$

(59)

Following, concerning the comparison with the current SoA, the total paid amount for a specific budget year ($t\in [1,I]$) is computed using equation 60:

$$R^{\left(t\right)}=\sum _{i+j=t}{y}_{i,j}\mathrm{with}t\in [1,I]$$

(60)

where the variable $y_{i,j}$ represents the single payments as indicated in Table 1 (see Section 2.1). Since the variables $y_{i,j}$ represent the observation values, they are not used as KPI to perform the comparison with the SoA. In the following Section, the comparison with the literature using the developed stochastic model for both residues will be shown, jointly with the obtained values of the $S\left(\mathit{\beta }\right)$ function for each residue, the dispersion parameters $\varphi$, the prediction of the claims reserves for each individual budget year $R^{\left(t\right)}$, and the entire claims reserve $R_{tot}$ to be set aside with the associated absolute errors.

3.1. Comparison with [22]

The data-set used in the work of [22] presents a Run-Off triangle filled in the context of “RC cars”, which is a well-known benchmark to test and validate the new methodology and theory of forecasting model especially in claims reserve estimation due to their small data dispersion and homogeneity. The authors applied the method previously described, using both residues defined in Section 2.4 to reproduce the results of [22]. The developed stochastic approach has been applied to predict the claims reserve $R^{\left(t\right)}$ for the generic budget year (t) using equation 58, and it has been compared with the one provided by the optimization of [22]. The results concerning the predicted claims reserve $R^{\left(t\right)}$ for the generic budget year (t), are depicted in the following two figures. Figure 1 shows the paid amount ($R^{\left(t\right)}$ in the equation 60) and the predicted claims reserve ($R^{\left(t\right)}$ in equation 58) relative to the comparison between [22] and our methodology using the residue defined in equation 46. Figure 2 shows the paid amount and the predicted claims reserve relative to the comparison between [22] and our methodology using the residue defined in equation 47. The obtained numerical results for the predicted claims reserve for a specific budget year are summarized in Table 3.

The obtained numerical results comprehending the value of the $S\left({\mathit{\beta }}^{*}\right)$ function, the corresponding dispersion parameter ($\varphi$), and the total predicted total claims reserve amount $R_{tot}$ are summarized in Table 4.

It is essential to mention that in both cases, i.e. in the estimation provided by the authors, the convergence of the GN algorithm is improved because the authors used the GA to calculate the initial guess of parameters $\mathit{\beta }$ (see Section 2.5). In contrast, concerning the work of [22], the algorithm’s convergence might converge to a local minimum without starting with an appropriate initial guess, leading to an increase in the errors associated with the $\mathit{\beta }$ parameters. Additionally, it is important to mention that, the computed values of the $S\left({\mathit{\beta }}^{*}\right)$ function, the dispersion parameter ($\varphi$), and the total predicted claims reserve amount $R_{tot}$ for the work of [22] have been obtained using the ${\mathit{\beta }}^{*}$ parameters resulted by their calculation. The errors on the predicted claims reserve for the specific year $R^{\left(t\right)}$ have been calculated using the error propagation by summing the errors relative to every single cell of the predicted values ($\Delta {\mu }_i$ - see equation 57). The error propagation theory has been applied to calculate the entire predicted claims reserve $R_{tot}$. Furthermore, using the two different definitions of the residues, the estimation of the claims reserve performed by the authors provides confidence intervals slightly lower than the ones provided by [22]. Even if the difference between the two residues is negligible for this case, the high dispersion and data aleatory in other contexts, such as healthcare, might produce more prominent confidence intervals.

3.2. Comparison with [23]

In this Section, the authors performed the validation approach of Section 3.1 with the Run-Off triangle presented in the work of [23]. In this case, [23] did not compute the errors associated with their stochastic model because they focused on developing a robustification of the CL method. However, the possibility to reproduce and improve the results presented in the work of [23] represents an excellent manner to validate further the methodology described in this paper since they present a low disperse data-set. Following the identical procedure of Section 3.1, the authors show the results of the developed stochastic model using the two definitions of the residues in equations 46 and 47, respectively. The results concerning the predicted claims reserve $R^{\left(t\right)}$ for the generic budget year (t), are depicted in the following two figures. Figure 3 shows the paid amount ($R^{\left(t\right)}$ in the equation 60) and the predicted claims reserve ($R^{\left(t\right)}$ in equation 58) relative to the comparison between [22] and our methodology using the residue defined in equation 46. Figure 4 shows the paid amount and the predicted claims reserve relative to the comparison between [23] and our methodology using the residue defined in equation 47. The obtained numerical results for the predicted claims reserve for a specific budget year are summarized in Table 5.

In this case, the authors provided the results concerning the value of the $S\left({\mathit{\beta }}^{*}\right)$ function, the corresponding dispersion parameter ($\varphi$), and the claims reserve amount applied to the Run-Off in the work of [23]. The results are shown in Table 6. Since [23] did not provide the errors associated with their stochastic model, the authors are unable to perform a consistent comparison with the work of literature, but they only provide the claims reserve amount coming from a robustification of the CL method.

Even in this case, both cases in the estimation provided by the authors, the convergence of the GN algorithm is improved because the authors used the GA to calculate the initial guess of parameters $\mathit{\beta }$ (see Section 2.5). Furthermore, using the two different definitions of the residues, the calculation performed by the authors provides confidence intervals for the claims reserve that include the values calculated by [23] using a completely different approach.

4. The health care study: the case of claims in the Tuscany Region

The main objective is to apply the developed methodology to predict the total claims reserve to be set aside by the insurance agency. Specifically, the authors worked on the health insurance context of claims concerning the Tuscany Region in Italy to provide results for the reserve amount required by the region to manage the yearly balance finance system in healthcare. It is essential to notice that the authors carried on the stochastic analysis on claims reserve employing the two definitions of residues described in Section 2.4, providing accordingly to the SoA a new methodology to deal with claims reserve and the corresponding errors. The novelty and the contribution of our research take credit not only from the improvements compared with the SoA but also in the healthcare context. In order to test our methodology in the healthcare insurance system of the Tuscany Region, the author proceeded in the following manner. Firstly, the authors performed a preliminary analysis of the available healthcare insurance data-set to understand the useful information that can be retrieved in order to compile the Run-Off triangle. Secondly, the authors compared the CL methodology with the developed stochastic approach. The CL is a deterministic method used in different contexts to compute the one-year reserve. Subsequently, the authors reproduced the deterministic result with the developed stochastic model. Finally, the authors found the optimal parameters of the ODP to perform prediction for the entire claims reserve to provide an outcome to the region, which could be used in the management, together with the corresponding confidence intervals.

4.1. Preliminary analysis

The authors worked with the Centro Gestione Rischio Clinico of the Tuscany Region database concerning the healthcare claims reported from 2010 and 2021. In this Section, the authors performed a preliminary analysis retrieving the following information:

• Time evolution of the paid claims;
• Single payment for years of the specific claim;
• Total paid amounts.

Starting from the provided dataset, it is possible to compile the following Run-Off triangle, which is the base of our research (Section 2.1). The rows of the Run-Off triangle are filled with the incurred years (i.e. from 2010 to 2021). The columns of the Run-Off triangle are filled with the progression of the payments towards the years. In this way, each cell represents the delayed payment with respect to the corresponding incurred year. Following the Run-Off triangle of the health insurance claims for the Tuscany Region is shown Table 7:

The variable $y_{i,j}$ in Table 7 identifies the paid amount of the payment for claims that have been incurred in the $i^{th}$ year with a delay occurred in the $j^{th}$ year. As previously discussed, authors present the deterministic analysis of the Run-Off triangle shown in Table 7 employing the well-known CL approach. This analysis represents the benchmark of our research, as it will be compared with the proposed stochastic method by the authors in the following Sections.

4.2. Chain Ladder Results

The CL method provides a mechanism to forecast future loss development trends based on historical patterns. It offers a deterministic perspective on reserving losses for the subsequent year. This Section outlines the results of the CL method applied to the previous Run-Off triangle. The approach has been comprehensively explained in the referenced studies [22,26]. The results obtained by applying the CL method to the health insurance of the Tuscany Region are summarised in Table 8.

In this case, the predicted claims reserve amount $R^{\left(t\right)}$ appearing in Table 8 has been calculated using equation 4, while the total predicted claims reserve amount for the Tuscany Region in the incurred years 2010-2021 is roughly: $R_{tot}=194.32\times {10}^6€$. Figure 5 shows the predicted claims reserve amount $R^{\left(t\right)}$ for every single incurred year (2010-2021) employing the CL methodology. Despite the simplicity of the CL method, several issues make it unable to perform our analysis. Indeed it cannot provide and predict changes in the model, and it can only predict the one-year estimate of the claims reserve. Therefore, a stochastic method is needed to provide more complete results for claims reserve predictions and confidence intervals. In the following Section 4.3, the authors present the developed stochastic approach with the corresponding numerical results compared to the CL approach.

4.3. Stochastic Approach Results

The results obtained by applying the developed stochastic approach to estimate the claims reserve for the Tuscany Region are presented in this Section. In this Section, the values and the corresponding errors of the parameters (${\mathit{\beta }}^{*}$) resulting from the developed stochastic approach are shown for the Tuscany Region insurance claims using both residue defined in equations 46 and 47. The procedure described in Section 2.4 and Section 2.5 has been applied to calculate the ${\mathit{\beta }}^{*}$ parameters needed to compile the Run-Off triangle. The GA has been applied firstly to glimpse a good initial guess of $\mathit{\beta }$ parameters, and subsequently, the GN algorithm has been applied to reach the stationary point ${\mathit{\beta }}^{*}$. The following Table 9 provides the results of the entire estimation for the $\mathit{\beta }$ parameters and the corresponding errors.

At this step, it is possible to compute the expectation value ${\mu }_{i,j}$ for each insurance year (2010-2021) and fulfil the Run-Off triangle as shown in Section 3. The following Table 10 and Table 11 represent the Run-Off triangles compiled with the estimation provided by the entire procedure described in Section 2.4 and Section 2.6.

Given the expectation values ${\mu }_{i,j}$ shown in the previous Table 10 and Table 11, the authors provided the following figures representing the paid amount and the predicted claims reserve $R^{\left(t\right)}$ for the generic budget year (t) from 2010 to 2031 in the Tuscany Region.

Referring to the previous Figure 6 and Figure 7, the Table 12 presents the complete set of values concerning the paid amount and the predicted claims reserve $R^{\left(t\right)}$ for the generic budget year (t) from 2010 to 2031 in the Tuscany Region.

In Table 13, the authors provided the values of the KPI resulting from the estimation of the total reserve to be set aside for the Tuscany Region insurance claims. Specifically, the value of the $S\left({\mathit{\beta }}^{*}\right)$ function, the corresponding dispersion parameter ($\varphi$), the predicted claims reserve amount for the 2021 incurred year $R^{\left(2021\right)}$, and the total predicted claims reserve $R_{tot}$ computed using equation 59.

The results in Table 13 illustrate the predictive accuracy of the stochastic model developed by the authors, which was previously validated as described in Section 3. It is important to note that the absolute errors on the predicted claims reserve are considerable off an interval of acceptability (i.e. the estimated errors are roughly 50 % of the predicted values for the claims reserve). The reason is due to the high data dispersion that the stochastic model is unable to overcome. Nevertheless, compared to the current SoA, the analysis conducted by the authors for the healthcare sector is entirely satisfactory, as it represents a notable initial approach to estimating the claims reserve in a highly heterogeneous and mutable context.

5. Conclusion

In this work, the authors contributed to creating a robust protocol to deal with claims reserve in the healthcare system context. Even if the SoA already presented a well-defined approach for claims reserve in the healthcare system, the wide variability, peculiarity, and aleatory made necessary reformulation and subsequent improvements of the stochastic methodology usually adopted. Therefore, the authors contributed to the SoA, improving the management of claims reserve in such a way that is easy to reproduce by the institutions which are the final users (i.e. the Tuscany Region) without increasing computational costs. Possible developments of the method concern the possibility of expanding the Run-Off tables, predicting the durability of the entire claims life, and defining a method to improve the error estimations also for highly dispersed databases such as the one of the healthcare context as treated in this work. Furthermore, it is necessary to create a unique protocol for data management to improve the quality and accessibility of the data and develop a robust analysis of the incurred velocity.