Heat Equation as a Tool for Outliers Mitigation in Run-Off Triangles for Valuing the Technical Provisions in Non-Life Insurance Business

: Estimating outstanding claims reserves in the non-life insurance business is often impaired by outlier-contaminated datasets. Widely used methods to eliminate outliers in non-life development triangles are either limiting the number of outliers by robust statistical methods or by change of development factors. However, the whole estimation process is likewise adversely affected so that: or (ii) estimation via the bootstrap method is ineligible. In this paper, the properties of the heat equation are examined to obtain an outlier smoothing technique for development triangles. The heat equation in two dimensions is being applied on an outlier contaminated dataset where no individual data are available. As a result, we introduce a new methodology to (i) treat outliers in non-life development triangles, (ii) keep the total sum of all triangle payments, and (iii) provide acceptable differences between the original and the backward estimated triangle. Consequently, the outlying values are eliminated and the resulting development triangle could be used as an input for any claims reserving method without a need for further robustiﬁcation or change of development factors. Additionally, the research on the application of heat equation in one dimension presented in this paper enables one to employ the bootstrap method using Pearson’s residuals in cases where the method was originally inapplicable due to development factors being lower than one.


Introduction
According to Solvency II regulations, the insurer must be able to estimate the future claims reserves as accurately as possible. The insurer who operates in the non-life insurance business often does not know the amount of the final claims for the year of the accident at the end of that year. It depends on the business line in the non-life insurance industry or the time duration of a claims settlement. Delays can occur due to the time lag between the occurrence of the accident and the appearance of the consequences of the event. Therefore, a run-off triangle can be considered to arrange the claims reserves. Most important is to estimate the outstanding claims reserve. Various methods can be used, the most popular one is the classic chain-ladder method (Verdonck et al. 2009).
"The chain ladder method is based on the assumption that the expectations underlying the columns and the rows in the run-off triangle are proportional" (Verdonck et al. 2009). Since the early 1990s, several articles have been published to incorporate the simple chain ladder method into the statistical framework, and consideration has been given to stochastic models that generate a chain ladder algorithm (Manolache 2019). Full stochastic models for the chain ladder method were published by (Liu and Verrall 2010;Mack 1993;Murphy 1994;Verdonck and Debruyne 2011;Verdonck et al. 2009) and other authors. Extended versions of these models have also been created. For example, (Peters et al. 2014) developed an extended class of model structures for the paid-incurred chain ladder models, where they developed exactly the Bayesian formulation of such models. (Wuthrich 2017) extended the chain-ladder method for claims reserving to include information about the properties of individual claims applying a neural network model. "The chain ladder method should only be used for large portfolios where consistency of the estimates is more important than unbiasedness and where all entries into the incurred loss triangle are rather reliable (no big relative chance fluctuations and/or errors)" (Bühlmann 2016, p. 7).
The Cape Cod method, which is also known as the Stanard-Buhlmann method, (Bühlmann and Straub 1983;Stanard 1985) was proposed to overcome some of the shortcomings of the chain ladder method (Saluz 2015). The a priori loss ratio in the Cape Cod method "is calculated as the weighted average of the chain ladder ultimate loss ratios across all years with the used premium as the weights" (Korn 2016, p. 1). Due to its simplicity and advantages over the chain ladder method, the Cape Cod method has become a proven method in practice (Saluz 2015). The Cape Cod method is a special case of the Generalized Cape Cod Methods addressed by (Gluck 1997;Korn 2016;Struzzieri et al. 1998).
A popular method that generates a simulated prediction distribution to obtain the standard errors of well-specified models is bootstrapping (Verdonck et al. 2009). This method has already been considered in the area of claims reserving by (England and Verrall 1999;Lowe 1994;Maciak et al. 2022;Zaçaj et al. 2022). Several authors ( Peremans et al. 2017;Verdonck et al. 2009) applied robust bootstrap procedures for the chain-ladder method.
One of the key decisions in estimating claims reserves is how to treat outliers. According to (Embrechts and Wüthrich 2022, p. 5) outliers in insurance typically are not data errors but large financial claims that are an important pricing component. (Verdonck and Van Wouwe 2011) proposed two techniques to detect and correct outliers in the bivariate chain-ladder method-the first technique was based on the bagplot to the bivariate dataset and the second one was the robust technique based on the MCD (Minimum Covariance Determinant). (Avanzi et al. 2022) extended their approach and also applied two alternative robust bivariate chain-ladder techniques to treat outlier-the first one was based on the outlyingness and the second technique was based on bagdistance, which is derived from the bagplot.
In relation to the above-mentioned, the objective of designing an in-house application (Barlak 2021) for computing non-life reserves using well-defined deterministic and stochastic methods (Avanzi et al. 2016;Badounas and Pitselis 2020;Brazauskas et al. 2009;Peremans et al. 2018;Verdonck and Van Wouwe 2011;Verdonck et al. 2009), and challenges associated with the lack of person-specific data, lead us to the design of a new method to treat outliers in non-life development triangles. By applying properties of the heat equation, outliers could be treated without changing the whole triangle payments total sum. Furthermore, a bootstrap method using residuals could be applied in some cases where it was originally impossible.
The heat equation is a partial differential equation that describes temperature changes in a given area over a period of time (Gorguis and Chan 2008). The one-dimensional heat equation was first studied by Fourier at the beginning of the 19th century (Cannon and Browder 1984). The heat equation has applications in various fields of science, one of the most important of them is the theory of heat conduction (Widder 1976). It has also been used for image enhancement (Black and Sapiro 1999;Buades et al. 2006) or for a detection of a pollution problem (El Badia and Ha-Duong 2002). (Itkin et al. 2021) applied multi-layer heat equations when solving financial problems and developing efficient algorithms for pricing barrier options for time-dependent one-factor short-rate models. At present, we are not aware of the use of the heat equation to treat outliers in non-life development triangles. In this paper, we fill this gap by proposing a new method to treat outliers, which is based on a heat equation.
Based on the above, we set the aim of the paper to design a new method to treat outliers in non-life development triangles.
The remainder of the paper is structured as follows. Section 2 outlines the theoretical basis of methods for the calculation of technical reserves. Selected Chain-ladder and Cape Cod deterministic methods (Section 2.1.1) with their stochastic adjustments (Section 2.1.2) and additional stochastic modification (Cowell 2009) (Section 2.1.3) are introduced in this section. Two different methods for treating outliers in 2-D (Sections 2.4.1 and 2.4.2) and an approach to the numerical solution of one-and two-dimensional heat equation (Sections 2.2 and 2.3) constitute the core part of this research. Finally, a method for adjusting development factors to be greater than one (>1) without changing the total sum in a triangle row is proposed in this section, with description and real-world examples being presented. Section 3 lists the results of the practical application of the heat equation when treating outliers in non-life development triangles. Section 5 summarizes the essential conclusions resulting from the research and presents the significant findings.

Methods for Technical Reserves Calculation
Hereby, we would like to recall some concepts that underlie the technical reserves calculation briefly. For a detailed description, the interested reader could refer to (Avanzi et al. 2016;Badounas and Pitselis 2020;Brazauskas et al. 2009;Peremans et al. 2018;Verdonck and Van Wouwe 2011;Verdonck et al. 2009).

Deterministic Methods
Let us start with a brief overview of the well-known deterministic methods used for the calculation of technical reserves in this paper. We use two deterministic methods for the mentioned purposes. Firstly, it is the Chain-ladder methodĈ CL i, for I as maximum number of years from an event of a claim and J stands for the total number of development years. Secondly, we have chosen the Cape-Cod method using the following formula for technical reserves computation The choice for a stochastic process was subjectively the simplest one, namely the bootstrap. We have employed two modifications of the deterministic methods. The first one uses residuals.

Bootstrap Method Using Residuals
In general, we have a triangle of cumulative payments (Table 1).
In the next step, we estimate development factors, which help us to calculate cumulative payments in the lower triangle (Table 2).
For the moment, we use the diagonal and apply the formulaĈ to estimate observed values of payments backwards (Table 3). We then get the upper triangle of both estimated and observed payments, which we will introduce into a bootstrap algorithm.  We compute the unscaled Pearson's residuals using r i,j = Using C i,j = ∑ j k=0 X i,k we easily get the formula from (Pesta 2011) where i + 1 ≤ I + 1 likewise and additionally C i,j−1 =Ĉ i,j−1 = 0 for j = 1. The next is the bootstrap algorithm itself, where for 1 ≤ b ≤ B (B stands for the total number of bootstrap cycles) following steps are performed:

1.
By using random sampling with replacement from the set of residuals from the upper triangle without the diagonal elements The new non-cumulative upper triangle is then computed using r i,j = 3. Non-cumulative upper triangle of the new "observed" payments is then used as an input to the classic deterministic method (in our case, the Chain-ladder or the Cape Cod) to get the vector of reserves (b) R.
After B simulations, we get B columns of reserves and their sums (Table 4).
Source: own construction.
We then easily get the n-th percentile from the sorted reserves sums.

Chain-Ladder Bootstrap Method Using Local Development Factors
The method proposed in (Cowell 2009). Let us denote the observed local development . At the same time, the first column in a triangle of cumulative payments is equal to the first column in the non-cumulative one. Therefore, C i,0 = X i,0 , ∀i ∈ {0, . . . , I}. We can write the cumulative triangle schematically as in Table 5.
Where each C i,j from the upper triangle can be expressed by the formulâ To compute the reserves column, we need to fill the local development factors into the lower triangle as well. The Chain-ladder method uses one estimated factor for each unoccupied cell in the column. In this case, the local development factor estimationλ i,j is obtained as a random sample with replacement from the set of all observed local development factors for a given development year {λ k,j }, where k ∈ {0, . . . , J − i}.
Let us have a look at an example (Table 6). −→ 285λ 2,3 − − → factors estimates will always consist of just one number as the set we are sampling is only a single-element one. The advantage of this approach is its simple application to the bootstrap. In every cycle b (for 1 ≤ b ≤ B), we fill the lower triangle with estimates of the local development factors. The reserves column is computed as Table 7). Table 7. Computed reserves for one bootstrap cycle applying the estimated local development factors.
Again, we get B vectors of reserves and the percentiles from their sorted sums.

Heat Equation in One Dimension
For the sake of simplicity, we take the plain form of the heat equation (Gurevich 2016) ∂u(x, t) ∂t where u(x, t) stands for a temperature in a point of one dimensional space x and in a specific point of time t for x ∈ [0, L]. The initial condition (so called Dirichlet boundary condition) will be the temperature of the whole interval [0, L] at t = 0 i.e., u( The boundary conditions represent the temperature of the interval boundaries (in this case points 0 and L) for the whole time period e.g., u(0, t) = u(L, t) = 0 ∀x > 0. Notice the heat diffusion in the Figure 1 (illustrated employing an online differential equation solver (Silvestre n.d.). Apparently, the initial condition is the whole function at time t = 0. In this case, we use the so-called Neumann boundary conditions ∂u(0,t) ∂x = ∂u(L,t) ∂x = 0, meaning that heat will neither leave nor enter the system on its boundaries. Notice the outlier in Figure 1a. After a very short period of time (in fact, almost instantly- Figure 1b), the heat equation diffuses high outlying temperature to the surroundings, while the function of temperature itself hardly changes. We aim to benefit from this property.

Numerical Solution of the Heat Equation in One Dimension
By means of dx the left side of the Formula (2) can be rewritten as ∂u ∂t Now, consider the continuous interval [0, L] as a discreet set of points {x i } n i=1 , where ∀i ∈ 1, . . . , n − 1 : x +1 − x i = dx. The right side of the Equation (2) can be equally rewritten as Putting (3) and (4) together we get a discreet form of (2) and, therefore The space discretization is illustrated in Table 8.
Source: own construction.
We can see temperatures in the discreet one-dimensional space for each time point in columns 1 to n. Columns 0 and n + 1 are boundary conditions. In this case, we do not prefer the heat to enter nor to leave the space; therefore, we simply set the column 0 equal to the column 1 (and, respectively, the column n + 1 equal to the column n).
When using the numerical solution of a partial differential equation, it is very important to regard the stability of the solution. In this case, the stability condition must be satisfied (Gurevich 2016) dt ≤ 1 2 dx 2 .

Heat Equation in Two Dimensions
Let us examine the heat equation in two dimensions where u(x, y, t) is temperature at point (x, y) in a two-dimensional space and time t. For our purpose, we can simplify our space to (x, y) ∈ Λ = [0, L] × [0, L]. Thus, the initial condition will be considered as function of temperature on the whole subspace Λ at t = 0, i.e., u(x, y, 0) = f (x, y) for x ∈ Λ. Boundary conditions will be represented by the temperature on the edges of the subset Λ at each time point. In the case of Dirichlet boundary conditions, we set u(0, y, t), u(L, y, t), u(x, 0, t), u(x, L, t) as constants. Respectively, using the Neumann boundary conditions as constants, we set the derivatives of the mentioned points, where again, temperature does not leave or enter the boundary of Λ.
The following Figure 2 shows escaping heat from a very thin plate in time.

Numerical Solution of the Heat Equation in Two Dimensions
Let us rewrite the equation in continuous form (8) to a discreet form as i.e., Discretization scheme in time t = τ is then (Table 9): White part of the Table 9 contains temperatures mapped from the set Λ in a discreet form. Initial condition represents temperature within the subspace Λ at τ = 0. Grey part of the table comprises boundary conditions. Finally, the temperature at time τ + 1 is computed by means of u τ i,j−1 , u τ i,j and u τ i,j+1 according to the Equation (10). Stability of the numerical solution will be conditioned by the form dt ≤ dx 2 dy 2 2(dx 2 + dy 2 ) .

Heat Equation Application
The heat equation has many applications, one of which is damaged painting restorations. For the majority of points in a painting, when taking a very small surrounding of a given point, the "difference" between it and other points is small, unlike some scratches (or naturally some edges of objects) where the "difference" is obviously bigger. Using the heat equation, we can smooth these differences at a cost of distorting the edges. Thus, the adjusted result can be blurred. Notice the woman's forehead in the Figure 3. In this case, the heat diffusion had a sufficient amount of time to soften the damage but not enough to make it too blurred. Small red squares on the Figure 3 mark the damaged areas. On the right side, we can see the effect of the heat equation application. This characteristic of the heat equation was a motivation to study its application also in non-life development triangles.

Practical Use-Triangle Transformation
Function of the non-cumulative payment development, depending on the development year, is generally decreasing. Thus, we cannot use the heat equation on a pure triangle. The first approach is to transform the non-cumulative payments triangle before it enters the heat equation. This way appears incorrect, but for the sake of research, we introduce this method, as well. In this case, from each payment, we subtract the median of its column. The resulting triangle is then the initial condition for the heat equation. For every whole row or column of a triangle in time τ in a form {ξ 1 , ξ 2 , . . . , ξ k−1 , ξ k } τ , a vector with boundary conditions will simply be {ξ 1 , ξ 1 , ξ 2 , . . . , ξ k−1 , ξ k , ξ k } τ , ensuring that no heat will enter or leave. (Nota bene: these boundary conditions are in Neumann's form). To secure stability, we set dt = 0.05 with dx = dy = 1 and the number of steps being at least two for letting the heat propagate to at least all eight adjacent cells (see Figure 4b).
(a) One step (b) Two steps Figure 4. Heat diffusion to adjacent points.
For our purposes, we set a maximum number of steps to four. After running the heat equation, it is important to transform the adjusted triangle back to the original form by adding medians to its values. At this point, we can run methods for reserves computation on a new adjusted triangle. In practice, however, an insurance company does not have a stable number of clients during accident years, for instance, which calls for a solution.

Practical Use-Heat Equation Adjustment
Instead of transforming the triangle, we try to transform the heat equation itself to be able to cope with systematic changes during the accident years, as well as the natural decline of payments during the development period. Let us rewrite the Formula (10) to We compare incomparable values in all four round brackets. It is similar to comparing today's 100 € to 100 € in 5 years. Therefore, we simply use the equivalent to the compounding and discounting approach.
Let us define horizontal factors between the non-cumulative values of payments aŝ h j = m X 0,j+1 X 0,j , . . . , X I−1,j+1 X I−1,j , X 0,j , . . . , X I−1,j , where function m(x, w) is a weighted median of values x j and of their corresponding weights w j . In case a local horizontal factor does not make sense (e.g., division by zero), the weighted median is computed without this local factor and its corresponding weight. We set the horizontal factor to 0.5 if it does not exist or is lower than zero. As a consequence, we preserve decreasing in the non-cumulative payments during the development period. Whether a constant is to be precisely determined in this case is up to the reader.
The vertical factors are set simply aŝ to eliminate possible edge problems during the accident years. X i,0 is the first possible non-cumulative payment (development year zero). Weighted medians are not used here because payments up to one year are the biggest carrier of information. In practice, it is not unusual that the outlier in the later development period clearly damages the ratio of payments between the accident years. The vertical factor is set to 1 if it does not exist or is lower than zero. As a result, we assume no systematic change between the following accident years. For both factors (horizontal and vertical), we set their distance from zero to be at least 0.001 to avoid unreasonably high multiplications. Then the adjusted step for the numerical solution is set as

Non-Positive Increments Estimates
The stochastic modification methods described in Section 2.1.2 are based on the computation of the residuals using the Formula (1). The square root in the denominator implies that the bootstrap method cannot be used when estimated increments are zero or lower (meaning that the development factors must be more than one).
For this case, we use the one-dimensional heat equation to try to adjust the development factors to values greater than one. The initial condition for the row of the triangle {r 1 , r 2 , . . . , r k−1 , r k } will be the row itself. The boundary condition in time τ will be {r 1 , r 1 , r 2 , . . . , r k−1 , r k , r k } τ . The time step is set to dt = 0.05 to satisfy stability condition (7), having dx = 1. The maximum number of steps is set to 8. If, after running the heat equation, any development factor is still one or lower than the bootstrap, the methods from Section 2.1.2 will not be used.

Practical Use-Triangle Transformation
Let us start with the triangle with an artificially created outlier, the number with bold in Table 10. The original non-outlier value was 15,888,572. After three steps 1 of the two-dimensional heat equation, we get the result displayed in Table 11 with reduced outlier value with bold. For the following results, it is important to be aware of the fact that the reserves were computed using the original vector of earned premium from (Gatialova 2010) for Tables 12 and 13.  34,130,722 CC 28,506,180 BCL 34,496,727 47,019,901 64,831,451 80,843,126 BCC 28,868,084 34,393,630 39,639,634 42,741,521 BCL_C 29,435,838 47,904,121 54,952,155 55,865,054 Source: own production.  31,616,200 CC 28,492,405 BCL 31,560,638 36,415,770 41,173,703 44,909,633 BCC 28,372,479 31,116,043 33,344,732 34,786,776 BCL_C 30,043,105 33,191,035 39,069,138 40,057,328 Source: own production.
As expected, besides the slightly modified point estimates, using the heat equation caused lower variance in the bootstrap methods.

Non-Positive Increments Estimates
Let us take the original triangle from (Gatialova 2010) (Table 14) and set one increment to −10,000. It is up to the reader why −500 is not enough in this case. After running the heat equation in one dimension (apparently on the first row of the triangle), we get the following adjusted triangle (Table 15). Whereas the sum of payments in each row before and after the adjustment is the same. Rounded development factors are summarized in Table 16.  (Table 14) 1.54711 1.06390 0.99995 1.00561 After (Table 15) 1.56507 1.07802 1.00202 1.00532 Source: own production.
As soon as all the development factors are greater than one, the adjusted triangle can be used as an input for the bootstrap method using residuals. For other methods, the original triangle was used. Afterwards, the reserves shown in Table 17 are obtained.

Practical Use-Heat Equation Adjustment
Let us introduce a real insurance company example (Motor insurance). As NBS is allowed to publish only aggregated data, the triangle is multiplied by a constant. Numbers are therefore different, but the ratios stay the same. For our purposes, it is also sufficient to show only a part of the triangle as shown in Table 18. In this case (systematic changes during accident years), it is not appropriate to use the triangle transformation technique. Some outliers can be smoothed; some others will appear. Figure 5 shows boxplots of the local development factors before and after the inappropriate adjustment. Using the heat equation adjustment technique from the Section 2.4.2, we get the boxplots of the local development factors as shown in Figure 5, where existing outliers are smoothed and the new ones do not appear again as they do in Figure 6. Both figures represent the boxplots of the local development factors for each development year. The Figure 5 shows how the development factors change after applying the heat equation on the original development triangle using an inappropriate triangle transformation technique. In Figure 6, the appropriate heat equation adjustment technique has been used.

Discussion
There are different approaches to mitigate outliers during the loss reserves estimation. The mitigation of outliers is usually dependent on a specific methodology (Bornhuetter and Ferguson 1972;Bühlmann and Straub 1983;Stanard 1985;Taylor 1977;Verdonck et al. 2009) for calculating reserves and the subsequent creation of a robust alternative.
Current research on the removal of outliers in contaminated datasets for claims reserves is quite limited as the traditional methods are preferably applied. In the light of the most recent works published in this topic (Avanzi et al. 2022b;Badounas et al. 2022), our method takes on a different approach. First, outlying values are detected and smoothed by the heat equation application. The resulting development triangle could then be used as an input for any method without the need for its robustification.
The future work shall address the appropriate setting for the parameters of the numerical solution of the heat equation, such as the step length or the maximum number of steps. The current parameters are set as constants, which seems to be sufficient for the demonstration of the methodology.

Conclusions
Outlier management in the non-life development triangles for calculation of technical reserves estimates has been thoroughly addressed in the existing literature. This article presents another way to approach the problem naturally occurring in the non-life development triangles. The heat equation and its properties introduced in the paper help address the lack of person-specific data while keeping the whole original triangle payments total sum and reasonable differences between the original and the backward estimated triangle. The method enables to smooth out the outliers across the neighboring years instead of their elimination which yields more realistic results. Hence, the unscaled Pearson's residuals for the stochastic modifications and the one-dimensional heat equation are used to adjust the development factors to values greater than one. Otherwise, the bootstrap method cannot be used. Introducing the one-dimension and the two-dimension heat equation offers an uncommon overview of the simplest stochastic process coupled with the physical equation and their potential in overcoming the challenges of the field.