Next Article in Journal
A Generalized Linear Model and Machine Learning Approach for Predicting the Frequency and Severity of Cargo Insurance in Thailand’s Border Trade Context
Next Article in Special Issue
Unified Spatial Clustering of Territory Risk to Uncover Impact of COVID-19 Pandemic on Major Coverages of Auto Insurance
Previous Article in Journal
Enhancing Sell-Type Home Reversion Products for Retirement Financing
Previous Article in Special Issue
Credibility Distribution Estimation with Weighted or Grouped Observations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Stochastic Claims Reserve in the Healthcare System: A Methodology Applied to Italian Data

1
Department of Management, School of Advanced Studies Sant’Anna, 56127 Pisa, Italy
2
The BioRobotics Institute, School of Advanced Studies Sant’Anna, 56127 Pisa, Italy
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Risks 2024, 12(2), 24; https://doi.org/10.3390/risks12020024
Submission received: 23 October 2023 / Revised: 24 January 2024 / Accepted: 26 January 2024 / Published: 29 January 2024

Abstract

:
One of the challenges in the healthcare sector is making accurate forecasts across insurance years for claims reserve. Healthcare claims present huge variability and heterogeneity influenced by random decisions of the courts and intrinsic characteristics of the damaged parties, which makes traditional methods for estimating reserves inadequate. We propose a new methodology to estimate claim reserves in the healthcare insurance system based on generalized linear models using the Overdispersed Poisson distribution function. In this context, we developed a method to estimate the parameters of the quasi-likelihood function using a Gauss–Newton algorithm optimized through a genetic algorithm. The genetic algorithm plays a crucial role in glimpsing the position of the global minimum to ensure a correct convergence of the Gauss–Newton method, where the choice of the initial guess is fundamental. This methodology is applied as a case study to the healthcare system of the Tuscany region. The results were validated by comparing them with state-of-the-art measurement of the confidence intervals of the Overdispersed Poisson distribution parameters with better outcomes. Hence, local healthcare authorities could use the proposed and improved methodology to allocate resources dedicated to healthcare and global management.

1. Introduction

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 USI (2000). This damage is closely correlated with civil liability in healthcare, which profoundly impacts the economy and society. Indeed, 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% OECD (2021). It should be noted that approximately 15% of hospital expenditure in high-income countries is used to correct complications of care and patient harm OECD (2018). However, in these countries which have efficient healthcare systems, 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 Mushinski et al. (2022)—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 Bertoli and Grembi (2018); Grembi and Garoupa (2013). When an adverse event occurs, Italian healthcare liability is imposed on healthcare organizations and, potentially, the healthcare professionals who committed the error. Public and private healthcare 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 that of employees. Other measures include the possibility for organizations to utilize cover funds from internal provisions to cover all or part of the health liability risk. The money to pay insurance premiums or to provide coverage funds is derived from the national health fund, which is then distributed to 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 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 functions 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 part, the authors illustrate an application of the developed methodology to the disposable data coming from the Tuscany region1 to show its capability compared to the available SoA; (4) finally, a 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 (R2022b release).

2. Structure and Main Ideas

The generalized linear model (GLM) has become a significant reality in claims reserve analysis in the last few years England and Verrall (1999, 2002); Kendall and Stuart (1967); Larsen (2007); Verrall (2000). 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 of 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 Run-Off triangles Schmidt (2006) 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) method Hindley (2017); Schmidt (2012); Verrall (1991). 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 England and Verrall (2002) 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 Hess and Schmidt (2002); Verdonck et al. (2009), and various stochastic claim reserving models have been developed for general insurance to extrapolate the tail factors.
The employment of stochastic models replicates the conventional reserve estimates produced by the CL methods, including parametric curves and smoothing techniques to reshape the Run-Off triangles, enabling 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 performance Renshaw and Verrall (1998). In this regard, the GLM has become, in the last years, a robust stochastic methodology for predicting claims reserve Hindley (2017); Nelder and Wedderburn (1972). The GLM represents a wide range of models that allow the definition and estimation of parameters by maximizing the likelihood function Nelder and Wedderburn (1972). 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 Wedderburn (1974, 1976), 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 Strascia and Tripodi (2018), the authors assume a logarithmic mean–variance relation, which provides the well-known Overdispersed Poisson (ODP) distribution model. The maximization of the ODP quasi-likelihood function is performed using the Gauss–Newton (GN) algorithm Fletcher (2000); Strascia and Tripodi (2018). 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 optimizing two different residuals. Consequently, it gives a better overview of the claims reserve in the healthcare sector, where the strong aleatory and peculiar nature of claims’ lifespan needs 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 on the data-sets present in the SoA, particularly from Strascia and Tripodi (2018) and Verdonck et al. (2009). These validation tests helped to understand whether the developed model outperforms the ones in the literature. Subsequently, the developed model was 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.

3. Theoretical Background

In this Section, the authors present the theoretical background for estimating claims reserve. Section 3.1 describes the CL approach and how to construct the Run-Off triangles that represent the applied basis of the developed methodology. Section 3.2 provides an introduction to GLM to understand the context in which our methodology is applied. Section 3.3 introduces the definition of the quasi-likelihood function to make predictions of the response variable defined in the GLM framework. Section 3.4 describes the implementation of the Gauss–Newton (GN) algorithm that maximizes the quasi-likelihood function and its limitation. Section 3.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, Section 3.6 provides a workflow for how to estimate the errors of the quasi-likelihood function.

3.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 reported 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, Run-Off triangles Schmidt (2006) 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 ith “incurred” year with i [ 1 , I ] , i.e., the year in which the claims happen, and the column corresponds to the jth “progression” year with j [ 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 ith year and effectively paid in the jth year. In the following, Table 1 shows a representation of a 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 Strascia and Tripodi (2018) and Schmidt (2006), the first step of the CL approach is the definition of a Run-Off triangle of the cumulative payments, which is shown in the following Equation (1):
C i , j = k = 0 j y i , k
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 C ^ i , J is obtained through the dilation factors f ^ j ( t ) (called link-ratio), according to Equation (2):
C ^ i , J ( t ) = C i , t i j = t i J 1 f ^ j ( t )
where the index t [ 1 , I ] indicates the generic budget year, while the term f ^ j ( t ) represents the link-ratio which is defined through the following Equation (3):
f ^ j ( t ) = k = 1 t j 1 C k , j + 1 k = 1 t j 1 C k , j
Hence, given Equations (2) and (3), the CL method is able to estimate the claims reserve R ^ i ( t ) given by Equation (4):
R ^ i ( t ) = C ^ i , J ( t ) C i , t 1
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.

3.2. Generalized Linear Model

The GLM was introduced in the work of Nelder and Wedderburn (1972) 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 some explicative variables or predictors with small variation ε with expectations value E [ ε ] = 0 and variance var [ ε ] = σ 2 . 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 relation between the response variable and the explicative variables is not linear;
  • The explicative variables Y do not take values from the interval ( , ) ;
  • The variance σ 2 is not constant;
  • It is not possible to assume the response variable Y following the normal distribution.
The GLM allows the user to overcome previous limitations by expanding the sets of applicability of the standard linear model and has became a common tool in statistical analysis and prediction theory. Indeed, the GLM needs 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 family 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 β as defined in Equation (5):
f ( y ; β , ϕ ) = e α ( ϕ ) y β g ( β ) + h ( y ) + α ( ϕ , y )
where y stands for a single realization of Y, while the variable ϕ represents the nuisance parameter2 and α ( ϕ ) is a positive definite function. In order to define the moments of the distributions, we introduce the likelihood L ( β | y ) and the log-likelihood L ( β | y ) functions as shown in Equation (6). In the following, the single realization (y) of the response variable (Y) refers to incremental payments as shown in Table 1.
L ( β | y ) = ln L ( β | y ) = ln i = 1 m f ( y i ; β , ϕ ) = i = 1 m f ( y i ; β , ϕ )
For the sake of simplicity, the authors moved from a matrix-like notation to a vector-like notation (with a single index). Therefore, the vectorial response variable y i , which is in one-to-one correspondence with the previous y i , j , is defined as in the following Table 2:
Based on the results of Kendall and Stuart (1967), it is possible to derive two crucial properties of the log-likelihood function shown in Equation (7):
E L β = 0 E 2 L β 2 = E L β 2
which provides the statement shown in Equation (8) for the generic exponential function (Equation (6)):
E L β = E α ( ϕ ) [ y g ( β ) ] = 0 E [ y ] μ = g ( β )
together with the quadratic expectation value:
E 2 L β 2 = E [ α ( ϕ ) g ( β ) ] E [ α ( ϕ ) g ( β ) ] = E [ α 2 ( ϕ ) ( y μ ) 2 ]
and hence,
g ( β ) = α ( ϕ ) E [ ( y μ ) 2 ]
which finally provides the expression for the variance of y (Equation (11)):
α ( ϕ ) var ( y ) = g ( β ) = V ( μ )
It is easy to prove that for a response variable Y normally distributed, the GLM reduces to the ordinary linear model.

3.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 (6)). In the work of Wedderburn (1974), a new methodology was proposed introducing the so-called quasi-likelihood function, which only needs to fix the relationship between the first two moments of the distribution, so, mean and variance. The mean–variance relation is entirely defined by the link function (g) as shown in Equation (12).
var ( y ) = 1 α ( ϕ ) V ( μ ) = 1 α ( ϕ ) g ( β ) = 1 α ( ϕ ) g ( ( g ) 1 ( μ ) )
In Equation (12), the term ( g ) 1 represents the inverse function of the first derivative of the link function (g) according to Equation (8).
Furthermore, we will assume, without any lack of generality, an exponential link function. Our proposed model robustly represents the well-known GLM with Overdispersed Poisson distribution (ODP). It is possible to prove that the overall function α ( ϕ ) 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 (12) can be written as indicated in Equation (13):
var ( y ) = ϕ g ( ( g ) 1 ( μ ) )
Consequently, for an exponential link function, the mean–variance relation is given by Equation (14):
E [ y ] = μ = 1 ϕ var ( y ) = V ( μ )
Considering that V ( μ ) = μ in Equation (14), the nuisance parameter is estimated by the Pearson statistics:
ϕ ϕ ^ = 1 m n i = 1 m ( y i μ i ) 2 V ( μ i ) ϕ = 1 m n i = 1 m ( y i μ i ) 2 μ i
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 , μ i ) is defined in Equation (16):
K ( y i , μ i ) μ i = y i μ i var ( y i ) = y i μ i ϕ V ( μ i )
or equivalently employing its integral form (Equation (17)):
K ( y i , μ i ) = y i μ i ϕ V ( μ i ) d μ i + c ( y i )
where the index i indicates the single observations. The variable c ( y i ) 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 Wedderburn (1974). Therefore, substituting Equation (14) into Equation (17) and integrating over the i-interval [ y i , μ i ] , we obtain the following expression (Equation (18)) for the quasi-likelihood function of the ODP model:
K ( y i , μ i ) = y i μ i y i μ i ϕ V ( μ i ) d μ i = y i μ i y i μ i ϕ μ i d μ i = 1 ϕ y i ln μ i y i μ i + y i
Finally, summing over the complete set of the Y realization, it is possible to obtain the complete quasi-likelihood function in terms of y, β , and ϕ :
K ( y , β , ϕ ) = 1 ϕ i = 1 m y i ln μ i y i μ i + y i
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 , μ 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 , μ 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 , μ i ) function. For clarity, the variable μ i will be indicated as μ i ( β ) because it represents the model that contains the parameters β that maximize the quasi-likelihood function. Additionally, the details of its convergence problem and a solution to this problem will be proposed for a better choice of initial guesses. In the following Sections, the authors will employ the vectorial notation for the β parameters for completeness.

3.4. Gauss–Newton Algorithm

Taking β = ( β 1 , , β n ) as the vector of parameters on which the response variable depends, to determine the β * vector that maximizes the quasi-likelihood function (see Equation (18)), we need to solve the following system of non-linear Equation (20):
K ( y i , μ i ( β ) ) β k = 0
An iterative algorithm is required to locate the kth parameters of the β vector, where k [ 1 , n ] , that solve the non-linear equation system shown in Equation (20). Therefore, the GN algorithm has been deployed based on the guidelines provided in the work of Wedderburn (1974). Given a set of generic functions (called residuals) { r i ( β ) } i = 1 m , the GN locates the set of β k parameters that minimize the function S ( β ) defined in Equation (21):
S ( β ) = i = 1 m r i 2 ( β )
where the specific choice of the residual function r i ( β ) is crucial in the entire process of parameter estimation. Furthermore, the authors tested two different r i ( β ) functions, which we will define later, providing different results. Starting from the quasi-likelihood function K ( y i , μ i ( β ) ) , the Taylor expansion around the parameters β is needed to minimize the function S ( β ) . Considering an expansion step δ given by Equation (22):
δ = β * β
The first-order Taylor expansion in δ can be written in Equation (23):
K ( y i , μ i ( β * ) ) = K ( y i , μ i ( β ) ) + K ( y i , μ i ( β ) ) β | β δ + o ( δ 2 )
It is important to note that, because the function K ( y i , μ i ( β ) ) is maximized by β * , its first derivative in Equation (23) is null when β = β * ; therefore, a maximization procedure must be developed to locate the δ * parameters. The GN works directly on the function S ( β ) defined in Equation (21), which has no triviality problem. Subsequently, it remains to define a consistent definition of the function r i ( β ) to start the optimization procedure properly and minimize the S ( β ) function through the GN methodology. In this Section, we propose two different definitions of residuals r i ( β ) , which allows us to work separately on forecasting the claims reserves 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 residuals, as characterised by Strascia and Tripodi (2018), concerns the direct optimization of the parameters inside the quasi-likelihood K ( y i , μ i ( β ) ) , and it is given by Equation (24):
r i ( β ) = y i ϕ K ( y i , μ i ( β ) ) r i ( β ) = 1 ϕ μ i ( β ) y i ln μ i ( β ) y i
The second formulation of the residuals, developed by the authors, concerns the optimization of the dispersion parameters ( ϕ -see Equation (15)), and it is given by the Equation (25):
r i ( β ) = y i μ i ( β ) μ i ( β )
Regarding the following steps, the authors will not use one of the previous definitions of residuals but will keep the discussion in general terms. The comparison of the two specific formulations of r i ( β ) (Equations (24) and (25)) will be discussed in Section 4 and Section 5, where the results related to comparison with the SoA, and the data from the Tuscany region are presented, respectively. The function S ( β ) does not present triviality problems compared to the quasi-likelihood function. Therefore, the first-order Taylor expansion of Equation (21) can be performed in δ (see Equation (22)) around the initial point β using a vectorial formalism which provides Equation (26):
S ( β * ) = S ( β ) + S ( β ) β | β δ + o ( δ 2 ) = S ( β ) + 2 r T ( β ) · r ( β ) β | β δ + o ( δ 2 )
The first derivative in Equation (26) represents the Jacobian matrix of the residuals computed in the initial point β . More precisely, the first derivative appearing in Equation (26) depends on the definition of residuals (see Equations (24) and (25)). Therefore, the single elements of the Jacobian matrix can be defined by the following Equations (27) and (28):
r i ( β ) β k | β = 1 ϕ μ i ( β ) β k | β y i μ i ( β ) μ i ( β ) J i k ( β ) ( see Equation ( 24 ) )
r i ( β ) β k | β = 1 2 μ i ( β ) β k | β y i + μ i ( β ) μ i ( β ) 3 2 J i k ( β ) ( see Equation ( 25 ) )
The term appearing on the left-hand side of Equations (27) and (28) represents the derivative of the ith-residual r i ( β ) space with respect to every β k -parameter. The variable J i k ( β ) appearing on the right-hand side of Equations (27) and (28) represents the i-row and k-column of the Jacobian matrix defined in Equation (26); thus, the Jacobian matrix has a dimension of m × n where m is the number of residuals and n represents the number of parameters (i.e., the length of the β vector). At this step, the previous Taylor expansion shown in Equation (26) can be rewritten as follows (Equation (29)):
S ( β * ) = S ( β ) + 2 r T ( β ) · J ( β ) · δ + o ( δ 2 ) = r T ( β ) · r ( β ) + 2 r T ( β ) · J ( β ) · δ + o ( δ 2 )
The term J ( β ) appearing in Equation (29) represents the Jacobian matrix in a vectorial formalism computed in the initial point β . The variable S ( β ) in Equation (29) is expressed by performing the dot product between the residuals r ( β ) and its vectorial transpose. For the sake of simplicity, the variables J ( β ) and r ( β ) will be indicated as J and r , respectively. At this point, it is possible to compute the gradient of the S ( β ) function in the stationary point β * as it appears in Equation (29):
S ( β ) β | β * = 2 J T r + 2 J T J δ + 2 r T J β | β δ
Since the gradient of S ( β ) is computed in the stationary point β * , the derivatives of this gradient are null. Additionally, Equation (30) can be further simplified by neglecting the term containing the derivative of the Jacobian matrix. Therefore, Equation (30) can be rewritten as follows:
S ( β ) β | β * = 2 J T r + 2 J T J δ = 0
Therefore, the δ can be recalculated as a function of β * (see Equation (22)) as follows:
δ = ( J T J ) 1 J T r β * = β ( J T J ) 1 J T r
In practice, the GN algorithm computes a β ( s ) value at every iterative step (s) starting from an initial guess of β ( 0 ) . Therefore, the δ ( s ) value is calculated using Equation (33):
δ ( s ) = ( J T ( β ( s ) ) · J ( β ( s ) ) ) 1 J T ( β ( s ) ) · r ( β ( s ) )
The iterative procedure described in Equation (33) is carried on until the following condition is verified (Equation (34)):
δ ( s ) β ( s + 1 ) β ( s ) < ϵ
where the variable ϵ represents an arbitrary user-defined tolerance which the GN algorithm must reach to stop the iterations. At this point, the GN algorithm has found the β parameters which minimize the function S ( β ) and therefore, the final value of β is given as follows (Equation (35)):
β * β ( s + 1 ) = β ( s ) ( J T ( β ( s ) ) · J ( β ( s ) ) ) 1 J T ( β ( s ) ) · r ( β ( s ) )
The described procedure represents an iterative methodology to locate the stationary point that maximizes the quasi-likelihood function K ( y i , μ i ( β ) ) , and it can be iterated until the δ reaches an arbitrary tolerance value ( ϵ = 10−15) defined by the user. For each step of the iterative method, the β is updated closer to the value that minimizes the S ( β ) function. The choice of a good initial guess for β is a crucial point for the iterative GN method to have success in locating the global minimum of the S ( β ) function. However, estimating a good initial guess of the multi-parametric function S ( β ) to force the convergence of the GN method to a global minimum requires additional concerns and tools that must be addressed.

3.4.1. The Hessian Modification

The convergence of the GN algorithm (Equation (35)) is not guaranteed for the residuals defined in Equations (24) and (25). The main reason that the GN algorithm fails to reach the convergence is due to the approximation introduced in Equation (31) where the derivative of the Jacobian matrix ( J ) with respect to β is ignored. In order to be able to apply Equation (35), the following approximation must hold:
J T J r T J β
In order to expect the convergence of the GN algorithm using Equation (35), the approximation highlighted in Equation (36) must be valid in one of the two cases: (1) the residuals r ( β ) are small compared to the derivatives of the Jacobian ( J / β ) around the stationary point β * ; (2) the residuals r ( β ) are blandly non-linear in such a way that the derivatives of the Jacobian are negligible Nocedal and Wright (2006). Therefore, in the case for which the approximation (36) cannot hold, the key equation of the GN algorithm (Equation (35)) must be modified. Starting from Equation (30), it can be rearranged in the following equation:
S ( β ) β | β * = 2 J T r + 2 J T J + 2 r T J β | β · δ S ( β * ) = S ( β ) + H S ( β ) · δ
The variable S ( β ) in Equation (37) represents the gradient of the function S ( β ) (see Equation (21)) computed in the initial point β ; the term H S ( β ) indicates the Hessian of the function S ( β ) computed in the initial point β . The Hessian is a matrix with a dimension of n × n where n represents the number of parameters (i.e., the length of the β vector), where each element of this matrix contains the second derivatives of the function S ( β ) . For simplicity, the gradient S ( β ) will be indicated as g , and the Hessian H S ( β ) will be indicated as H . Since the gradient S ( β * ) is computed in the stationary point β * , the derivatives of this gradient are null. Therefore, Equation (37) can be rewritten as follows:
S ( β * ) = g + H · δ = 0 δ = H 1 · g
Therefore, the δ can be recalculated as function of β * (see Equation (22)) as follows:
β * = β H 1 · g
In this case, the GN algorithm computes a β ( s ) value at every iterative step (s) starting from an initial guess of β ( 0 ) . Therefore, the δ ( s ) value is calculated using the following Equation (40):
δ ( s ) = H ( β ( s ) ) 1 · g ( β ( s ) )
The iterative procedure described in Equation (40) is carried on until the condition δ < ϵ is verified. The variable ϵ represents an arbitrary user-defined tolerance which the GN algorithm must reach to stop the iterations. At this point, the GN algorithm has found the β parameters which minimize the function S ( β ) ; therefore, the final value of β is given as follows (35):
β * β ( s + 1 ) = β ( s ) H ( β ( s ) ) 1 · g ( β ( s ) )
The described procedure represents an alternative iterative methodology to locate the stationary point ( β * ) that maximizes the quasi-likelihood function K ( y i , μ i ( β ) ) , and it can be adopted when the previous methodology (see Section 3.4) fails completely to locate the stationary point.

3.4.2. Limits of the GN Algorithm

The convergence of the GN algorithm using the Jacobian matrix (see Section 3.4) and its Hessian modification (see Section 3.4.1) is not always guaranteed, as proved and discussed by Nocedal and Wright (2006) and Mascarenhas (2013). Even local convergence could fail for the following reasons: (1) the initial guess β ( 0 ) is far from the stationary point β * ; (2) the J T J matrix described in Equation (32) is ill conditioned; (3) the Hessian matrix H in Equation (39) is ill conditioned. Therefore, in order to perform claims reserve estimation in this framework, it is necessary to address and solve the problem of GN convergence and provide a method to ensure efficiency in finding the global minimum by the algorithm. Local convergence allows us to predict acceptable results for the claims reserve but could provide huge errors associated with the estimated β parameters. Thus, providing an instrument that guarantees the GN algorithm’s global convergence is necessary to offer good results for claims reserve estimation and associated errors without increasing computational costs. In this regard, the author aims to derive a good initial guess of the β parameters to prevent divergence or local convergence of the GN algorithm. The genetic algorithm (GA) has been employed to help locate a good initial guess of the β parameters to solve the convergence problem.

3.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 β 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: Conn et al. (1991, 1997); Goldberg (1989). The algorithm, at each iteration, creates a sequence of new populations of individuals starting from the current generation. In order to avoid a huge computational cost, the author decided to implement a population of 2000 individuals, with 20,000 generations as the maximum number of generations, and 5000 generations as the maximum stall generation number for which the fitness function remains constant during the optimization.
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 residuals defined in Equation (24) or Equation (25). 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 β parameters by combining them as an entry of a pair of parents-crossover or a single parent-mutation. The parent-crossover procedure follows the well-known technique described by Haupt and Haupt (2004), where two elite individuals are selected to generate the children for the subsequent generation. Let us suppose that the two parents are labelled as a and b, where their corresponding β -parameters are defined as the following vectors [ β 1 a , β 2 a ] and [ β 1 b , β 2 b ], respectively. The procedure involves the random selection of a single β i -parameter which will serve as the crossover point, and subsequently, a random number ( η ), which ranges from 0 to 1, is introduced to determine the new β i -parameters for the children in the next generation. The calculation of the new β i -parameter is given by the following Equation (42):
β i 1 = ( 1 η ) β i a + η β i b
β i 2 = ( 1 η ) β i b + η β i a
The other β -parameter, which is β 2 in our example, is directly inherited from each respective parent chromosome. At this point, suppose that the two children are labelled as c and d; the complete set of β -values for the next generation is given by the following vectors: [ β i 1 c , β 2 c ] and [ β i 2 d , β 2 d ]. The described methodology combines the characteristics of both parent chromosomes using a single crossover point jointly with a η random value to obtain a new set of chromosomes. In contrast, the parent-mutation involves a random change of a single β i -parameter within a user-defined range [ β i m i n , β i m a x ]. The new parent parameter for the next generation is given by the following relationship (44):
β i = β i m i n + η ( β i m a x β i m i n ) with η [ 0 , 1 ]
Additionally, the GA excludes individuals with low fitness who are not participating in the selection process. Finally, the stopping criterion of the GA is given by reaching a certain deviation threshold (fixed at 10 10 ) within the β -parameters in the fitness values for all the selected population Conn et al. (1991, 1997). The employment of the GA enforces the stability of the GN algorithm by increasing the convergence speed with a consequent reduction of the number of iterations needed to null the gradient represented in Equation (38) and ensuring that the global minimum is reached within the user-defined domain of the β parameter.

3.6. Expectation Values and Error Estimation

Following the approach described in Section 3.4, the parameterization of the β vector for the quasi-likelihood function defined in Equation (17) must be carried out. The parameterization described by Strascia and Tripodi (2018) has been followed to define the β vector as shown in Equation (45):
β = ( c , a 1 , , a n , b 1 , , b m )
where the authors have used a different definition of the indexes n and m with respect to Section 3.4, which is the one that will be used in the following. Therefore, the index is k [ 1 , n + m + 1 ] and it labels the k-parameters of the β vector. The variable c represents the intercept, the variable a n i (with n i [ 1 , n ] ) is associated with the n i t h rows of the Run-Off triangles, i.e., the incurred year, and the variable b m i (with m i [ 1 , m ] ) corresponds to the m i t h columns of the Run-Off triangles, i.e., the progression year. In this case, this notation has the advantage that the subscript index (i) indicates the ith residual, 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 β is represented by Equation (46) and indicates the expectation value of the y i observable:
E [ y i ] = μ i ( β ) = e c if n i = m i = 1 e c + a n i if m i = 1 and n i 1 e c + b m i if n i = 1 and m i 1 e c + a n i + b m i otherwise
Once the ODP framework has been chosen, the methodology discussed in Section 3.4 can be applied. The main goal is to find the best estimate of the β parameters associated with the generic μ i model (Equation (46)) 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 maximize the quasi-likelihood function (17) by substituting it in Equation (46). In this manner, it is possible to obtain Equation (47):
K ( y i , μ i ( β ) ) = 1 ϕ y i ( c + a n i + b m i ) y i log y i e c + a n i + b m i + y i
At this step, the function S ( β ) defined in Equation (21) has been built with the two residuals r i ( β ) (Equations (24) and (25)) defined in Section 3.4 and it depends directly on the parameter β . In the ODP framework, the residuals are defined as the following Equations (48) and (49):
r i ( β ) = 1 ϕ e c + a n i + b m i + y i log y i y i ( c + a n i + b m i ) ( Equation ( 24 ) )
r i ( β ) = y i e c + a n i + b m i e 1 2 ( c + a n i + b m i ) ( Equation ( 25 ) )
At this step, considering the total number of residuals (N) which corresponds to the number of observations, the Jacobian matrix ( J ) to apply the GN algorithm (see Equation (35)) is defined as follows (Equation (50)):
J = r 1 c r 1 a 1 r 1 a n r 1 b 1 r 1 b m r i c r i a 1 r i a n r i b 1 r i b m r N c r N a 1 r N a n r N b 1 r N b m
Starting from the initial guess β ( 0 ) provided by the GA (Section 3.5), the GN algorithm updates the Jacobian matrix (Equation (50)) using one of the two residuals defined in Equations (48) or (49) to compute iteratively the parameters β * which minimize the function S ( β ) through Equation (35). During the iterative process using the Jacobian matrix, the relation (36) must hold to reach the convergence of the calculus. If Equation (36) cannot hold, the Hessian modification of the GN algorithm (Section 3.4.1) must be applied. In this case, the gradient ( g ) and the Hessian ( H ) of the function S ( β ) are defined as follows (for simplicity, the variable S ( β ) is denoted as S):
g = S c , S a 1 , , S a n , S b 1 , , S b m T
H = 2 S c 2 2 S c a 1 2 S c a n 2 S c b 1 2 S c b m 2 S a 1 c 2 S a 1 2 2 S a 1 a n 2 S a 1 b 1 2 S a 1 b m 2 S a n c 2 S a n a 1 2 S a n 2 2 S a n b 1 2 S a n b m 2 S b 1 c 2 S b 1 a 1 2 S b 1 a n 2 S b 1 2 2 S b 1 b m 2 S b m c 2 S b m a 1 2 S b m a n 2 S b m b 1 2 S b m 2
At that point, it is possible to fill the entire Run-Off table with the complete set of estimations given by Equation (53):
μ i ( β * ) = e c * if n i = m i = 1 e c * + a n i * if m i = 1 and n i 1 e c * + b m i * if n i = 1 and m i 1 e c * + a n i * + b m i * otherwise
Concurrently with the estimation of the expectation values, the corresponding errors on the β * parameters have been computed starting from the Hessian matrix H ( K ( y , β * , ϕ ) ) of the complete quasi-likelihood function defined in Equation (19), where each element of the Hessian matrix ( H j k ) is computed using Equation (54):
H j k = 2 K ( y , β * , ϕ ) β j β k
Therefore, the variance ( δ β 2 ) associated with each β * parameter is computed by performing the inverse of the Hessian matrix ( H ), and then taking the diagonal elements of the matrix changing the sign as shown in Equation (55). For the sake of simplicity, the β * term is indicated as β , and the variable H indicates the Hessian of the complete quasi-likelihood function H ( K ( y , β * , ϕ ) ) .
δ β 2 = δ c 2 , δ a 1 2 , , δ a n 2 , δ b 1 2 , , δ b m 2 = diag ( H 1 )
In order to provide the complete estimation of the standard errors for the expectation values μ i ( β ) , the authors applied the theory of error propagation. Therefore, they computed the differential of the generic expectation values μ i ( β ) (see Equation (53)), which is given by Equation (56):
d μ i = d ( e c + a n i + b m i ) = e c + a n i + b m i ( d c + d a n i + d b m i )
Hence, substituting the parameter errors δ c , δ a n i , and δ b m i with the corresponding differentials in Equation (56):
δ μ i = e c + a n i + b m i ( δ c + δ a n i + δ b m i )
The variance δ μ i of the expectation values is given by the squared sum of the single term of the previous Equation (57):
δ μ i 2 = e 2 ( c + a n i + b m i ) ( δ c 2 + δ a n i 2 + δ b m i 2 + 2 δ c δ a n i + 2 δ c δ b m i + 2 δ a n i δ b m i )
where the mixed terms are negligible for independent variables. In this framework, the standard errors of the expectation value Δ μ i are defined as follows (59):
Δ μ i = e c + a n i + b m i d c 2 + d a n i 2 + d b m i 2
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 Section 3.3, Section 3.4 and Section 3.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 3.2, composed by expectation value and corresponding standard error. The authors are now able to present, in the next Section, a validation of the developed approach by comparison with the available SoA.

4. 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 Strascia and Tripodi (2018), which primarily inspired the authors because it presents a similar methodology, presenting a sweetened data-set from a context different from healthcare (i.e., RC cars in Italy). The data-set of Strascia and Tripodi (2018) provides a complete calculation of the expectation values with the corresponding errors to the estimated parameters ( β ). The second comparison concerns the work of Verdonck et al. (2009), 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 3.4, Section 3.5 and Section 3.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 residuals defined, respectively, in Equations (48) and (49). The key performance indicators (KPIs) used to compare our methodology with the ones proposed in the literature are the function S ( β ) computed as shown in Equation (21) using the values of the β vector that minimize the residuals after applying the GN algorithm ( β * ), the dispersion parameter ( ϕ ) computed as shown in Equation (15), and the claims reserve amount ( R t o t ) that must be set aside including the corresponding errors of the expectation value.
The developed stochastic model allows us to compute claims reserve differently with respect to the CL method described in Section 3.1. After computing the β * parameters of the stationary point through the GN algorithm, the expectation values shown in Equation (53) 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 3.
The coloured cells in “light blue” as shown in Table 3 represent the expectation values μ i , j with i + j I coming from the application of the GN algorithm. In other words, they correspond to the generic model μ i ( β * ) indicated in Section 3.6, for which they have been used to construct the residuals r i ( β ) in the Equations (48) and (49). Additionally, these expectation values ( μ i , j ) are important to understand the dispersion of the observable data ( y i , j in Table 1) and therefore compute the dispersion parameter ( ϕ ). In contrast, the coloured cells in “light orange” as shown in Table 3 represent the expectation values μ i , j with i + j > I that are estimated using Equation (53), and they are being used to predict the claims reserve for the generic budget year ( t > I ) through Equation (60):
R ( t ) = i + j = t μ i , j with t [ I + 1 , I + J ]
In order to predict the entire claims reserve R t o t , i.e., the total amount to be set aside, it is necessary to sum over every budget year t [ I + 1 , I + J ] , through Equation (61):
R t o t = t = I + 1 I + J R ( t )
Furthermore, concerning the comparison with the current SoA, the total paid amount for a specific budget year ( t [ 1 , I ] ) is computed using Equation (62):
R ( t ) = i + j = t y i , j with t [ 1 , I ]
where the variable y i , j represents the single payments as indicated in Table 1 (see Section 3.1). Since the variables y i , j represent the observation values, they are not used as KPIs to perform the comparison with the SoA. In the following Section, the comparison with the literature using the developed stochastic model for both residuals will be shown, jointly with the obtained values of the S ( β ) function for each residual, the dispersion parameters ϕ , the prediction of the claims reserves for each individual budget year R ( t ) , and the entire claims reserve R t o t to be set aside with the associated absolute errors. For this comparison, the authors used the data-sets of Strascia and Tripodi (2018) and Verdonck et al. (2009), both of which are accessible in their entirety within their respective articles. For the sake of simplicity, in the following the authors will only report the results of the comparison between the results in the aforementioned works and those given by their methodology, described in Section 3. The source data, arranged in Run-Off triangles, were selected for comparison as they relate to incremental payments over time in the claims reserve context, though not appropriate to the healthcare context specifically addressed in this paper. The comparison, based on the specific KPIs described above, is made to show how with the methodology proposed by the authors the results of the SoA can be reproduced with greater absolute and percentage accuracy, thanks to the implementation of the GA which is fundamental in the choice of the initial guess, a step not envisaged in the proposed literature.

4.1. Comparison with Strascia and Tripodi (2018)

The data-set used in the work of Strascia and Tripodi (2018) 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 residuals defined in Section 3.4, to reproduce the results of Strascia and Tripodi (2018). The developed stochastic approach has been applied to predict the claims reserve R ( t ) for the generic budget year (t) using Equation (60), and it has been compared with the one provided by the optimization of Strascia and Tripodi (2018). The results concerning the predicted claims reserve R ( t ) for the generic budget year (t), are depicted in the following two figures. Figure 1 shows the paid amount ( R ( t ) in Equation (62)) and the predicted claims reserve ( R ( t ) in Equation (60)) relative to the comparison between Strascia and Tripodi (2018) and our methodology using the residual defined in Equation (48). Figure 2 shows the paid amount and the predicted claims reserve relative to the comparison between Strascia and Tripodi (2018) and our methodology using the residual defined in Equation (49). The obtained numerical results for the predicted claims reserve for a specific budget year are summarized in Table 4.
The obtained numerical results describing the value of the S ( β * ) function, the corresponding dispersion parameter ( ϕ ), and the total predicted total claims reserve amount R t o t are summarized in Table 5.
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 β (see Section 3.5). In contrast, concerning the work of Strascia and Tripodi (2018), 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 β parameters. Additionally, it is important to mention that the computed values of the S ( β * ) function, the dispersion parameter ( ϕ ), and the total predicted claims reserve amount R t o t for the work of Strascia and Tripodi (2018) have been obtained using the β * parameters resulting from their calculation. The errors in the predicted claims reserve for the specific year R ( t ) have been calculated using the error propagation by summing the errors relative to every single cell of the predicted values ( Δ μ i —see Equation (59)). The error propagation theory has been applied to calculate the entire predicted claims reserve R t o t . Furthermore, using the two different definitions of the residuals, the estimation of the claims reserve performed by the authors provides confidence intervals slightly lower than the ones provided by Strascia and Tripodi (2018). Even if the difference between the two residuals is negligible for this case, the high dispersion and aleatory data in other contexts, such as healthcare, might produce more prominent confidence intervals.

4.2. Comparison with Verdonck et al. (2009)

In this Section, the authors performed the validation approach of Section 4.1 with the Run-Off triangle presented in the work of Verdonck et al. (2009). In this case, Verdonck et al. (2009) 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 Verdonck et al. (2009) represents an excellent manner to validate further the methodology described in this paper since they present a low-dispersion data-set. Following an identical procedure to the one of Section 4.1, the authors show the results of the developed stochastic model using the two definitions of the residuals in Equations (48) and (49), respectively. The results concerning the predicted claims reserve R ( t ) for the generic budget year (t), are depicted in the following two figures. Figure 3 shows the paid amount ( R ( t ) in Equation (62)) and the predicted claims reserve ( R ( t ) in Equation (60)) relative to the comparison between Strascia and Tripodi (2018) and our methodology using the residual defined in Equation (48). Figure 4 shows the paid amount and the predicted claims reserve relative to the comparison between Verdonck et al. (2009) and our methodology using the residual defined in Equation (49). The obtained numerical results for the predicted claims reserve for a specific budget year are summarized in Table 6.
In this case, the authors provided the results concerning the value of the S ( β * ) function, the corresponding dispersion parameter ( ϕ ), and the claims reserve amount applied to the Run-Off in the work of Verdonck et al. (2009). The results are shown in Table 7. Since Verdonck et al. (2009) did not provide the errors associated with their stochastic model, the authors are unable to perform a consistent comparison with this work in the literature, but they only provide the claims reserve amount coming from a robustification of the CL method.
Even in this case, and in 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 β (see Section 3.5). Furthermore, using the two different definitions of the residuals, the calculation performed by the authors provides confidence intervals for the claims reserve that include the values calculated by Verdonck et al. (2009) using a completely different approach.

5. The Healthcare 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 an 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 out a stochastic analysis on claims reserve employing the two definitions of residuals described in Section 3.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 authors 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 could 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 a 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 for the region which could be used in management, together with the corresponding confidence intervals.

5.1. Preliminary Analysis

The authors worked with the Centro Gestione Rischio Clinico of the Tuscany region database concerning the healthcare claims reported from 2010 to 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 data-set, it is possible to compile the following Run-Off triangle, which is the base of our research (Section 3.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. The Run-Off triangle of health insurance claims for the Tuscany region is shown in Table 8:
The variable y i , j in Table 8 identifies the paid amount of the payment for claims that have been incurred in the ith year with a delay occurring in the jth year. As previously discussed, the authors present the deterministic analysis of the Run-Off triangle shown in Table 8 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.

5.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. This approach has been comprehensively explained in the referenced studies Strascia and Tripodi (2018); Verrall (1991). The results obtained by applying the CL method to the health insurance context of the Tuscany region are summarized in Table 9.
In this case, the predicted claims reserve amount R ( t ) appearing in Table 9 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 t o t = 194.32 × 106 €. Figure 5 shows the predicted claims reserve amount R ( t ) for every single incurred year (2010–2021) employing the CL methodology. Despite the simplicity of the CL method, several issues make it unsuitable for performing 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 5.3, the authors present the developed stochastic approach with the corresponding numerical results compared to the CL approach.

5.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. Here, the values and the corresponding errors of the parameters ( β * ) resulting from the developed stochastic approach are shown for the Tuscany region insurance claims using both residual defined in Equations (48) and (49). The procedure described in Section 3.4 and Section 3.5 has been applied to calculate the β * parameters needed to compile the Run-Off triangle. The GA has been applied firstly to obtain a good initial guess of β parameters, and subsequently, the GN algorithm has been applied to reach the stationary point β * . The following Table 10 provides the results of the entire estimation for the β parameters and the corresponding errors.
At this step, it is possible to compute the expectation value μ i , j for each insurance year (2010–2021) and fulfil the Run-Off triangle as shown in Section 4. The following Table 11 and Table 12 represent the Run-Off triangles compiled with the estimation provided by the entire procedure described in Section 3.4 and Section 3.6.
Given the expectation values μ i , j shown in the previous Table 11 and Table 12, the authors provide the following figures representing the paid amount and the predicted claims reserve R ( t ) for the generic budget year (t) from 2010 to 2031 in the Tuscany region.
Referring to the previous Figure 6 and Figure 7, Table 13 presents the complete set of values concerning the paid amount and the predicted claims reserve R ( t ) for the generic budget year (t) from 2010 to 2031 in the Tuscany region.
In Table 14, the authors provided the values of the KPIs resulting from the estimation of the total reserve to be set aside for the Tuscany region insurance claims. Specifically, the value of the S ( β * ) function, the corresponding dispersion parameter ( ϕ ), the predicted claims reserve amount for the 2021 incurred year R ( 2021 ) , and the total predicted claims reserve R t o t are computed using Equation (61).
The results in Table 14 illustrate the predictive accuracy of the stochastic model developed by the authors, which was previously validated as described in Section 4. It is important to note that the absolute errors on the predicted claims reserve are considerably off an interval of acceptability (i.e., the estimated errors are roughly 50% of the predicted values for the claims reserve). The reason for this is 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.

6. Conclusions

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 nature of claims made reformulation and subsequent improvements of the stochastic methodology usually adopted necessary. 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 that are the final users (i.e., the Tuscany region) without increasing computational costs. This study wants to prove how it is possible to handle healthcare claims using the well-known GLM with an Overdispersed Poisson function with an improvement in the estimation of the parameters of the quasi-likelihood function, with an initial guess optimized by a genetic algorithm which glimpses the position of the global minimum solving the convergence of the Gauss–Newton algorithm. This step is crucial for the peculiarity and variability of healthcare databases. Possible developments of the study concern the possibility of expanding the Run-Off tables, predicting the durability of the entire claim life, and defining a method to better improve error estimations also for highly dispersed databases such as the one treated in this work in the healthcare context. 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.
[custom] References

Author Contributions

Conceptualization, G.C. and M.V.; Methodology, C.M. and A.D.; Software, C.M. and A.D.; Formal analysis, C.M. and A.D.; Writing—review and editing, C.M., A.D. and A.V.; Supervision, G.C. and M.V. All authors have read and agreed to the published version of the manuscript.

Funding

The research project developed in this paper is part of a broader research work funded by the Tuscany Region.

Data Availability Statement

Data were obtained from Centro Gestione Rischio Clinico of the Tuscany region and are unavailable due to privacy restrictions.

Acknowledgments

The authors sincerely thank the Tuscany region and, particularly, Federico Gelli and Michela Maielli; the Centro Gestione Rischio Clinico, and in particular Elisa Flore, Francesca Guerrini, and Pasquale Macrì; Nino Savelli, and Gian Paolo Clemente for the important feedback and support, and Chiara Seghieri from MeS laboratory of Sant’Anna School for Advanced Studies.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this article.

Notes

1
Data in the present paper are smoothed for policy and legal motivations.
2
A nuisance parameter is any unspecified parameter necessary to ensure that the model describes the system adequately. In our case, it represents the dispersion of the measured data.

References

  1. Bertoli, Paola, and Veronica Grembi. 2018. Courts, scheduled damages, and medical malpractice insurance. Empirical Economics 55: 831–54. [Google Scholar] [CrossRef]
  2. Conn, Andrew, Nicholas Gould, and Philippe Toint. 1991. A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds. SIAM Journal on Numerical Analysis 28: 2. [Google Scholar] [CrossRef]
  3. Conn, Andrew, Nicholas Gould, and Philippe Toint. 1997. A globally convergent augmented lagrangian barrier algorithm for optimization with general inequality constraints and simple bounds. Mathematics of Computation 66: 261–68. [Google Scholar] [CrossRef]
  4. England, Peter, and Richard Verrall. 1999. Analytic and bootstrap estimates of prediction errors in claims reserving. Insurance: Mathematics and Economics 25: 281–93. [Google Scholar] [CrossRef]
  5. England, Peter D., and Richard J. Verrall. 2002. Stochastic claims reserving in general insurance. British Actuarial Journal 8: 443–544. [Google Scholar] [CrossRef]
  6. Fletcher, Roger. 2000. Practical Methods of Optimization, 2nd ed. London: Wiley. [Google Scholar] [CrossRef]
  7. Goldberg, David E. 1989. Genetic Algorithms in Search, Optimization and Machine Learning. Boston: Addison-Wesley Longman Publishing Co., Inc. [Google Scholar]
  8. Grembi, Veronica, and Nuno Garoupa. 2013. Delays in medical malpractice litigation in civil law jurisdictions: Some evidence from the italian court of cassation. Health Economics, Policy and Law 8: 423–52. [Google Scholar] [CrossRef] [PubMed]
  9. Haupt, Randy L., and Sue Ellen Haupt. 2004. Practical Genetic Algorithms. London: John Wiley & Sons. [Google Scholar] [CrossRef]
  10. Hess, Klaus Th., and Klaus D. Schmidt. 2002. A comparison of models for the chain–ladder method. Insurance: Mathematics and Economics 31: 351–64. [Google Scholar] [CrossRef]
  11. Hindley, David. 2017. Claims Reserving in General Insurance. Cambridge: Cambridge University Press. [Google Scholar] [CrossRef]
  12. Institute of Medicine (U.S.) Committee on Quality of Health Care in America. 2000. To Err Is Human-Building a Safer Health System. Washington, DC: National Academies Press. [Google Scholar]
  13. Kendall, Maurice George, and Alan Stuart. 1967. The Advanced Theory of Statistics, 2nd ed. New York: Hafner. [Google Scholar]
  14. Larsen, Christian Roholte. 2007. An individual claims reserving model. ASTIN Bulletin 37: 113–32. [Google Scholar] [CrossRef]
  15. Mascarenhas, Walter. 2013. The divergence of the bfgs and gauss-newton methods. Mathematical Programming 147: 1–2. [Google Scholar] [CrossRef]
  16. Mushinski, David, Sammy Zahran, and Aanston Frazier. 2022. Physician behaviour, malpractice risk and defensive medicine: An investigation of cesarean deliveries. Health Economics, Policy and Law 17: 247–65. [Google Scholar] [CrossRef] [PubMed]
  17. Nelder, John Ashworth, and Robert William Maclagan Wedderburn. 1972. Generalized Linear Model. Journal of the Royal Statistical Society: Series A 135: 370. [Google Scholar] [CrossRef]
  18. Nocedal, Jorge, and Stephen J. Wright. 2006. Numerical Optimization, 2nd ed. New York: Springer. [Google Scholar]
  19. OECD. 2018. Delivering Quality Health Services: A Global Imperative for Universal Health Coverage. Paris: OECD Publishing. [Google Scholar]
  20. OECD. 2021. State of Health in the EU: Italy. OECD, European Commission & European Observatory on Health Systems and Policies. Paris: OECD Publishing. [Google Scholar]
  21. Renshaw, Arthur E., and Richard J. Verrall. 1998. A stochastic model underlying the chain-ladder technique. British Actuarial Journal 4: 903–23. [Google Scholar] [CrossRef]
  22. Schmidt, Klaus D. 2006. Methods and models of loss reserving based on run-off triangles: A unifying survey. Casualty Actuarial Society Forum 2006: 269–317. [Google Scholar]
  23. Schmidt, Klaus D. 2012. Loss prediction based on run-off triangles. AStA Advances in Statistical Analysis 96: 265–310. [Google Scholar] [CrossRef]
  24. Strascia, Stefano, and Agostino Tripodi. 2018. Overdispersed-poisson model in claims reserving: Closed tool for one-year volatility in glm framework. Risks 6: 139. [Google Scholar] [CrossRef]
  25. Verdonck, Tim, Martine Van Wouwe, and Jan Dhaene. 2009. A robustification of the chain-ladder method. North American Actuarial Journal 13: 280–98. [Google Scholar] [CrossRef]
  26. Verrall, Richard J. 1991. Chain ladder and maximum likelihood. Journal of the Institute of Actuaries (1886–1994) 118: 489–99. [Google Scholar] [CrossRef]
  27. Verrall, Richard J. 2000. An investigation into stochastic claims reserving models and the chain-ladder technique. Insurance: Mathematics and Economics 26: 91–99. [Google Scholar] [CrossRef]
  28. Wedderburn, Robert W. M. 1974. Quasi-likelihood functions, generalized linear models, and the gauss-newton method. Biometrika 61: 439–47. [Google Scholar] [CrossRef]
  29. Wedderburn, Robert W. M. 1976. On the existence and uniqueness of the maximum likelihood estimates for certain generalized linear models. Biometrika 63: 27–32. [Google Scholar] [CrossRef]
Figure 1. Computed reserve amount between our approach and the one provided by Strascia and Tripodi (2018), using the residual defined in Equation (48) for the GN algorithm. The “blue” bars indicate the paid amounts computed using Equation (62), while the “yellow” and the “orange” error-bars indicate the predicted claims reserve computed using Equation (60) for Strascia and Tripodi (2018) and our work, respectively.
Figure 1. Computed reserve amount between our approach and the one provided by Strascia and Tripodi (2018), using the residual defined in Equation (48) for the GN algorithm. The “blue” bars indicate the paid amounts computed using Equation (62), while the “yellow” and the “orange” error-bars indicate the predicted claims reserve computed using Equation (60) for Strascia and Tripodi (2018) and our work, respectively.
Risks 12 00024 g001
Figure 2. Computed reserve amount between our approach and the one provided by Strascia and Tripodi (2018), using the residual defined in Equation (49) for the GN algorithm. The “blue” bars indicate the paid amounts computed using Equation (62), while the “yellow” and the “orange” error-bars indicate the predicted claims reserve computed using Equation (60) for Strascia and Tripodi (2018) and our work, respectively.
Figure 2. Computed reserve amount between our approach and the one provided by Strascia and Tripodi (2018), using the residual defined in Equation (49) for the GN algorithm. The “blue” bars indicate the paid amounts computed using Equation (62), while the “yellow” and the “orange” error-bars indicate the predicted claims reserve computed using Equation (60) for Strascia and Tripodi (2018) and our work, respectively.
Risks 12 00024 g002
Figure 3. Computed reserve amount between our approach and the one provided by Verdonck et al. (2009), using the residual defined in Equation (48) for the GN algorithm. The “blue” bars indicate the paid amounts computed using Equation (62), while the “yellow” and the “orange” error-bars indicate the predicted claims reserve computed using Equation (60) for Verdonck et al. (2009) and our work, respectively.
Figure 3. Computed reserve amount between our approach and the one provided by Verdonck et al. (2009), using the residual defined in Equation (48) for the GN algorithm. The “blue” bars indicate the paid amounts computed using Equation (62), while the “yellow” and the “orange” error-bars indicate the predicted claims reserve computed using Equation (60) for Verdonck et al. (2009) and our work, respectively.
Risks 12 00024 g003
Figure 4. Computed reserve amount between our approach and the one provided by Verdonck et al. (2009), using the residual defined in Equation (49) for the GN algorithm. The “blue” bars indicate the paid amounts computed using Equation (62), while the “yellow” and the “orange” error-bars indicate the predicted claims reserve computed using Equation (60) for Verdonck et al. (2009) and our work, respectively.
Figure 4. Computed reserve amount between our approach and the one provided by Verdonck et al. (2009), using the residual defined in Equation (49) for the GN algorithm. The “blue” bars indicate the paid amounts computed using Equation (62), while the “yellow” and the “orange” error-bars indicate the predicted claims reserve computed using Equation (60) for Verdonck et al. (2009) and our work, respectively.
Risks 12 00024 g004
Figure 5. Predicted claims reserve amount R ( t ) for every single incurred year (2010–2021) using the CL method.
Figure 5. Predicted claims reserve amount R ( t ) for every single incurred year (2010–2021) using the CL method.
Risks 12 00024 g005
Figure 6. Paid amount and the predicted claims reserve R ( t ) for the generic budget year (t) from 2010 to 2031 in the Tuscany region computed using the residual defined in Equation (48).
Figure 6. Paid amount and the predicted claims reserve R ( t ) for the generic budget year (t) from 2010 to 2031 in the Tuscany region computed using the residual defined in Equation (48).
Risks 12 00024 g006
Figure 7. Paid amount and the predicted claims reserve R ( t ) for the generic budget year (t) from 2010 to 2031 in the Tuscany region computed using the residual defined in Equation (49).
Figure 7. Paid amount and the predicted claims reserve R ( t ) for the generic budget year (t) from 2010 to 2031 in the Tuscany region computed using the residual defined in Equation (49).
Risks 12 00024 g007
Table 1. Example of Run-Off triangle. The variable y i , j indicates the paid amount relative to a claim that has occurred in the ith incurred year where the effective payment has been delayed in the jth year.
Table 1. Example of Run-Off triangle. The variable y i , j indicates the paid amount relative to a claim that has occurred in the ith incurred year where the effective payment has been delayed in the jth year.
i / j 01j J 1 J
1 y 1 , 0 y 1 , 1 y 1 , j y 1 , J 1 y 1 , J
2 y 2 , 0 y 2 , 1 y 2 , j y 2 , J 1
⋮.
i y i , 0 y i , 1 y i , j
I 1 y I 1 , 0 y I 1 , 1
I y I , 0
Table 2. One-to-one correspondence between matrix-like and vector-like notation.
Table 2. One-to-one correspondence between matrix-like and vector-like notation.
Response VariableEntries
Matrix-like Notation
(Table 1)
y i , j y 1 , 1 , y 1 , 2 , , y k , l , , y m , n
Vector-like Notation y i y 1 , y 2 , , y k , , y m
Table 3. Complete set of expectation values μ i , j computed using the developed stochastic model.
Table 3. Complete set of expectation values μ i , j computed using the developed stochastic model.
i / j 01j J 1 J
1 μ 1 , 0 μ 1 , 1 μ 1 , j μ 1 , J 1 μ 1 , J
2 μ 2 , 0 μ 2 , 1 μ 2 , j μ 2 , J 1 μ 2 , J
i μ i , 0 μ i , 1 μ i , j μ i , J 1 μ i , J
⋮.
I 1 μ I 1 , 0 μ I 1 , 1 μ I 1 , j μ I 1 , J 1 μ I 1 , J
I μ I , 0 μ I , 1 μ I , j μ I , J 1 μ I , J
Table 4. Comparison between our developed methodology and the work of Strascia and Tripodi (2018) in predicting claims reserve R ( t ) for the specific budget year t [ I + 1 , I + J ] for the two residuals.
Table 4. Comparison between our developed methodology and the work of Strascia and Tripodi (2018) in predicting claims reserve R ( t ) for the specific budget year t [ I + 1 , I + J ] for the two residuals.
Budget
Year (t)
Paid ( 10 3 €)
Literature
R ( t ) ( 10 3 €)
Literature
R ( t ) ( 10 3 €)
Residual (Equation (48))
R ( t ) ( 10 3 €)
Residual (Equation (49))
022.60-22.60 ± 1.4322.60 ± 1.43
162.32-62.32 ± 5.7062.32 ± 5.65
2101.93-101.93 ± 9.25101.93 ± 9.03
3124.59-124.59 ± 12.12124.59 ± 11.57
4152.04-152.04 ± 14.82152.04 ± 13.85
5188.65-188.65 ± 17.33188.65 ± 16.03
6185.31-185.31 ± 19.72185.31 ± 18.45
7203.38-203.38 ± 21.77203.38 ± 20.40
8213.67-213.67 ± 23.03213.67 ± 21.37
9207.65-207.65 ± 23.39207.65 ± 21.63
10197.67-197.67 ± 23.64197.67 ± 21.84
11184.68-184.68 ± 24.99184.68 ± 23.10
12194.08-194.08 ± 29.10194.08 ± 27.00
13-177.71 ± 41.39183.87 ± 27.96179.53 ± 25.78
14-139.04 ± 33.41146.17 ± 22.90140.69 ± 20.80
15-112.39 ± 27.91119.75 ± 19.38113.99 ± 17.43
16-93.69 ± 24.06101.15 ± 16.9895.37 ± 15.13
17-80.55 ± 21.3685.66 ± 14.8981.93 ± 13.48
18-66.73 ± 18.3271.52 ± 12.9067.97 ± 11.61
19-52.05 ± 14.8457.26 ± 10.7253.14 ± 9.43
20-38.71 ± 11.5042.30 ± 8.2139.44 ± 7.26
21-29.33 ± 9.0532.14 ± 6.4429.97 ± 5.71
22-23.89 ± 7.6226.14 ± 5.4124.38 ± 4.78
23-18.78 ± 6.3920.36 ± 4.5019.19 ± 4.00
24-12.95 ± 4.9313.69 ± 3.4012.85 ± 3.01
Table 5. Comparison between the work of Strascia and Tripodi (2018) and our proposed methodology concerning the S ( β * ) function (Equation (21)), the dispersion parameter ϕ (Equation (15)), and the total predicted claims reserve amount R t o t (Equation (61)) for both residuals (Section 3.4).
Table 5. Comparison between the work of Strascia and Tripodi (2018) and our proposed methodology concerning the S ( β * ) function (Equation (21)), the dispersion parameter ϕ (Equation (15)), and the total predicted claims reserve amount R t o t (Equation (61)) for both residuals (Section 3.4).
Residual
(Equation (48))
Reference S ( β * ) ϕ R t o t ( 10 5 €)
Strascia and Tripodi (2018) 3.65 × 10 5 4.11 × 10 2 8.46 ± 2.20
Our work 3.00 × 10 5 4.53 × 10 2 9.00 ± 1.54
Residual
(Equation (49))
Reference S ( β * ) ϕ R t o t ( 10 5 €)
Strascia and Tripodi (2018) 2.71 × 10 4 4.11 × 10 2 8.46 ± 2.20
Our work 2.68 × 10 4 4.06 × 10 2 8.58 ± 1.38
Table 6. Comparison between our developed methodology and the work of Verdonck et al. (2009) in predicting claims reserve R ( t ) for the specific budget year t [ I + 1 , I + J ] for the two residuals.
Table 6. Comparison between our developed methodology and the work of Verdonck et al. (2009) in predicting claims reserve R ( t ) for the specific budget year t [ I + 1 , I + J ] for the two residuals.
Budget
Year (t)
Paid ( 10 6 €)
Literature
R ( t ) ( 10 6 €)
Literature
R ( t ) ( 10 6 €)
Residual (Equation (48))
R ( t ) ( 10 6 €)
Residual (Equation (49))
0135.34-135.34 ± 3.05135.34 ± 2.92
1216.03-216.03 ± 7.56216.03 ± 7.38
2294.31-294.31 ± 11.12294.31 ± 10.84
3353.38-353.38 ± 14.22353.38 ± 13.84
4431.01-431.01 ± 17.14431.01 ± 16.69
5463.80-463.80 ± 19.19463.80 ± 18.65
6478.03-478.03 ± 20.99478.03 ± 20.33
7493.10-493.10 ± 22.89493.10 ± 22.20
8507.59-507.59 ± 25.35507.59 ± 24.53
9535.26-535.26 ± 28.75535.26 ± 27.84
10-401.96403.74 ± 23.07402.54 ± 22.35
11-309.07310.82 ± 18.74309.57 ± 18.13
12-236.98238.43 ± 15.25237.42 ± 14.75
13-178.16179.39 ± 12.26178.56 ± 11.86
14-129.86131.02 ± 9.57130.23 ± 9.24
15-92.7893.62 ± 7.4593.08 ± 7.19
16-62.8563.36 ± 5.6163.08 ± 5.43
17-36.6337.03 ± 3.8236.83 ± 3.69
18-15.1214.87 ± 1.9315.11 ± 1.90
Table 7. Comparison between the work of Verdonck et al. (2009) and our proposed methodology concerning the S ( β * ) function (Equation (21)), the dispersion parameter ϕ (Equation (15)), and the total predicted claims reserve amount R t o t (Equation (61)) for both residuals (Section 3.4).
Table 7. Comparison between the work of Verdonck et al. (2009) and our proposed methodology concerning the S ( β * ) function (Equation (21)), the dispersion parameter ϕ (Equation (15)), and the total predicted claims reserve amount R t o t (Equation (61)) for both residuals (Section 3.4).
Residual
(Equation (48))
Reference S ( β * ) ϕ R t o t ( 10 9 €)
Verdonck et al. (2009)--1.46
Our work 8.57 × 10 6 2.02 × 10 5 1.47 ± 0.09
Residual
(Equation (49))
Reference S ( β * ) ϕ R t o t ( 10 9 €)
Verdonck et al. (2009)--1.46
Our work 6.81 × 10 6 1.89 × 10 5 1.47 ± 0.09
Table 8. Paid amount of the Run-Off triangle for health insurance claims for the Tuscany region from 2010 to 2021 (106 €).
Table 8. Paid amount of the Run-Off triangle for health insurance claims for the Tuscany region from 2010 to 2021 (106 €).
y i , j 01234567891011
20102.805.9410.15.004.064.282.741.982.641.061.201.27
20112.1311.36.151.792.824.295.267.764.913.351.43
20120.957.617.193.254.773.085.035.653.641.74
20130.728.665.681.633.444.064.893.952.50
20140.979.767.676.494.075.664.413.90
20151.5011.47.173.956.525.767.87
20161.155.149.149.743.332.95
20170.196.968.925.773.63
20180.105.909.363.68
20192.355.5510.9
20200.276.57
20210.81
Table 9. Predicted claims reserve amount R ( t ) for every single incurred year ( t ) using the CL method.
Table 9. Predicted claims reserve amount R ( t ) for every single incurred year ( t ) using the CL method.
t R ( t ) (€)
2011 1.56 × 10 6
2012 2.59 × 10 6
2013 3.97 × 10 6
2014 9.19 × 10 6
2015 16.83 × 10 6
2016 19.17 × 10 6
2017 22.24 × 10 6
2018 23.21 × 10 6
2019 34.07 × 10 6
2020 29.78 × 10 6
2021 31.70 × 10 6
Table 10. Estimated parameters ( β * ) using both residuals defined in Equations (48) and (49).
Table 10. Estimated parameters ( β * ) using both residuals defined in Equations (48) and (49).
Residual (Equation (48))Residual (Equation  (49))
VariableValueErrorVariableValueError
c14.1830.221c14.1090.226
a 1 0.1750.155 a 1 0.1970.156
a 2 −0.0220.164 a 2 0.0020.165
a 3 −0.1030.170 a 3 −0.1140.172
a 4 0.1260.165 a 4 0.1400.166
a 5 0.2860.165 a 5 0.3130.166
a 6 0.1910.177 a 6 0.1780.179
a 7 0.0480.193 a 7 0.0550.195
a 8 −0.0110.209 a 8 −0.0470.213
a 9 0.1640.217 a 9 0.1900.216
a 10 −0.1580.313 a 10 −0.1910.320
a 11 −0.5720.889 a 11 −0.4980.883
b 1 1.6550.205 b 1 1.7010.210
b 2 1.6820.207 b 2 1.7400.212
b 3 1.2440.220 b 3 1.2470.226
b 4 0.9940.233 b 4 1.0340.237
b 5 1.0260.237 b 5 1.0750.241
b 6 1.2160.238 b 6 1.2480.242
b 7 1.2390.248 b 7 1.2520.253
b 8 0.8940.285 b 8 0.9250.289
b 9 0.4300.355 b 9 0.4180.364
b 10 −0.1860.522 b 10 −0.1240.519
b 11 −0.1250.724 b 11 −0.0510.719
Table 11. Expectation values μ i , j (106 €) computed using the residual defined in Equation (48).
Table 11. Expectation values μ i , j (106 €) computed using the residual defined in Equation (48).
μ i , j 01234567891011
20101.447.557.765.013.904.034.874.983.532.221.201.27
20111.729.009.245.974.654.805.805.944.202.641.431.52
20121.417.397.594.903.823.944.774.873.452.171.171.25
20131.306.817.004.523.523.634.394.493.182.001.081.15
20141.648.578.805.684.424.575.525.654.002.521.361.44
20151.9210.110.36.675.195.366.486.634.702.951.601.70
20161.759.149.396.064.724.885.896.034.272.691.451.54
20171.527.938.145.264.094.235.115.233.702.331.261.34
20181.437.477.684.963.863.994.824.933.492.201.191.26
20191.708.909.145.904.604.755.745.874.162.621.411.50
20201.236.456.634.283.333.444.164.253.011.901.021.09
20210.814.264.382.832.202.272.752.811.991.250.680.72
Table 12. Expectation values μ i , j (106 €) computed using the residual defined in Equation (49).
Table 12. Expectation values μ i , j (106 €) computed using the residual defined in Equation (49).
μ i , j 01234567891011
20101.347.357.644.673.773.934.674.693.382.041.181.27
20111.638.959.305.684.594.795.695.714.122.481.441.55
20121.347.367.664.683.783.944.684.703.392.041.191.28
20131.206.566.824.173.373.514.174.193.021.821.061.14
20141.548.458.795.374.344.525.375.403.892.341.361.47
20151.8310.010.46.385.165.376.386.414.622.791.621.74
20161.608.789.135.584.514.705.585.614.042.431.421.52
20171.427.778.074.933.994.154.934.963.582.151.251.35
20181.287.017.294.453.603.754.464.483.231.941.131.22
20191.628.889.245.644.564.755.645.674.092.461.431.54
20201.116.076.313.863.123.253.863.882.791.680.981.05
20210.814.474.642.842.292.392.842.852.061.240.720.77
Table 13. Representation of the paid amount and predicted claims reserve R ( t ) for the generic budget year provided with the two different definitions of the residuals. From the budget year 2022 to 2031, the values of R ( t ) are fully predicted.
Table 13. Representation of the paid amount and predicted claims reserve R ( t ) for the generic budget year provided with the two different definitions of the residuals. From the budget year 2022 to 2031, the values of R ( t ) are fully predicted.
Budget Year (t) R ( t ) (106 €)
Paid Amount
R ( t ) (106 €)
Residual (Equation (48))
R ( t ) (106 €)
Residual (Equation (49))
20102.801.44 ± 0.321.34 ± 0.30
20118.089.27 ± 2.748.98 ± 2.71
201222.318.2 ± 5.7917.9 ± 5.83
201319.522.9 ± 7.6122.5 ± 7.63
201422.725.9 ± 8.7625.2 ± 8.69
201527.331.1 ± 10.630.3 ± 10.5
201633.638.6 ± 13.337.8 ± 13.3
201732.844.9 ± 15.743.7 ± 15.5
201843.747.7 ± 17.046.2 ± 16.8
201955.648.8 ± 17.947.0 ± 17.5
202046.650.5 ± 19.148.7 ± 18.8
202147.350.5 ± 20.748.8 ± 20.3
2022-46.1 ± 21.244.6 ± 21.0
2023-37.9 ± 18.136.7 ± 17.9
2024-30.7 ± 14.629.4 ± 14.3
2025-25.7 ± 12.624.8 ± 12.4
2026-21.3 ± 11.120.5 ± 10.9
2027-16.2 ± 9.2715.5 ± 9.08
2028-11.0 ± 7.1010.6 ± 6.97
2029-6.56 ± 4.756.39 ± 4.69
2030-3.78 ± 3.083.76 ± 3.06
2031-1.76 ± 1.601.77 ± 1.62
Table 14. Results concerning the Tuscany region analysis providing the values of the KPIs such as S ( β * ) function (Equation (21)), the dispersion parameter ϕ (Equation (15)), the predicted claims reserve amount for the 2021 incurred year R ( 2021 ) , and the total predicted claims reserve R t o t .
Table 14. Results concerning the Tuscany region analysis providing the values of the KPIs such as S ( β * ) function (Equation (21)), the dispersion parameter ϕ (Equation (15)), the predicted claims reserve amount for the 2021 incurred year R ( 2021 ) , and the total predicted claims reserve R t o t .
S ( β * ) ϕ R ( 2021 ) (106 €) R tot (106 €)
Measured--47.3-
CL (Table 9)--31.7-
Residual (Equation (48))6.78 × 1036.05 × 10550.5 ± 20.7201 ± 103
Residual (Equation (49))3.27 × 1075.94 × 10548.8 ± 20.3194 ± 102
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Mazzi, C.; Damone, A.; Vandelli, A.; Ciuti, G.; Vainieri, M. Stochastic Claims Reserve in the Healthcare System: A Methodology Applied to Italian Data. Risks 2024, 12, 24. https://doi.org/10.3390/risks12020024

AMA Style

Mazzi C, Damone A, Vandelli A, Ciuti G, Vainieri M. Stochastic Claims Reserve in the Healthcare System: A Methodology Applied to Italian Data. Risks. 2024; 12(2):24. https://doi.org/10.3390/risks12020024

Chicago/Turabian Style

Mazzi, Claudio, Angelo Damone, Andrea Vandelli, Gastone Ciuti, and Milena Vainieri. 2024. "Stochastic Claims Reserve in the Healthcare System: A Methodology Applied to Italian Data" Risks 12, no. 2: 24. https://doi.org/10.3390/risks12020024

APA Style

Mazzi, C., Damone, A., Vandelli, A., Ciuti, G., & Vainieri, M. (2024). Stochastic Claims Reserve in the Healthcare System: A Methodology Applied to Italian Data. Risks, 12(2), 24. https://doi.org/10.3390/risks12020024

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop