An individual claims history simulation machine

Gabrielli, Andrea; Wüthrich, Mario V. Article An individual claims history simulation machine Risks Provided in Cooperation with: MDPI – Multidisciplinary Digital Publishing Institute, Basel Suggested Citation: Gabrielli, Andrea; Wüthrich, Mario V. (2018) : An individual claims history simulation machine, Risks, ISSN 2227-9091, MDPI, Basel, Vol. 6, Iss. 2, pp. 1-32, https://doi.org/10.3390/risks6020029 This Version is available at: https://hdl.handle.net/10419/195821 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ risks Article An Individual Claims History Simulation Machine Andrea Gabrielli * and Mario V. Wüthrich RiskLab, Department of Mathematics, ETH Zürich, Rämistrasse 101, 8092 Zürich, Switzerland; [email protected] *Correspondence: andr[email protected] Received: 5 March 2018; Accepted: 27 March 2018; Published: 30 March 2018   Abstract: The aim of this project is to develop a stochastic simulation machine that generates individual claims histories of non-life insurance claims. This simulation machine is based on neural networks to incorporate individual claims feature information. We provide a fully calibrated stochastic scenario generator that is based on real non-life insurance data. This stochastic simulation machine allows everyone to simulate their own synthetic insurance portfolio of individual claims histories and back-test thier preferred claims reserving method. Keywords: claims reserving; individual claims; claims cash flows; micro-level stochastic reserving; loss reserving; claims simulation; neural network reserving; individual claims features; individual claims covariates; chain-ladder 1. Introduction The aim of this project is to develop a stochastic simulation machine that generates individual claims histories of non-life insurance claims. These individual claims histories should depend on individual claims feature information such as the line of business concerned, the claims code involved or the age of the injured. This feature information should influence the reporting delay of the individual claim, the claim amount paid, its individual cash flow pattern as well as its settlement delay. The resulting (simulated) individual claims histories should be as ‘realistic’ as possible so that they may reflect a real insurance claims portfolio. These simulated claims then allow us to back-test classical aggregate claims reserving methods—such as the chain-ladder method—as well as to develop new claims reserving methods which are based on individual claims histories. The latter has become increasingly popular in actuarial science, see Antonio and Plat (2014), Hiabu et al. (2016), Jessen et al. (2011), Martínez-Miranda et al. (2015), Pigeon et al. (2013), Taylor et al. (2008), Verrall and Wüthrich (2016) and Wüthrich (2018a) for recent developments. A main shortcoming in this field of research is that there is no publicly available individual claims history data. Therefore, there is no possibility to back-test the proposed individual claims reserving methods. For this reason, we believe that this project is very beneficial to the actuarial community because it provides a common ground and publicly available (synthetic) data for research in the field of individual claims reserving. This paper is divided into four sections. In this first section we describe the general idea of the simulation machine as well as the chosen data used for model calibration. In Section 2we describe the design of our individual claims history simulation machine using neural networks. Section 3focuses on the calibration of these neural networks. In Section 4we carry out a use test by comparing the real data to the synthetically generated data in a chain-ladder claims reserving analysis. Appendix A presents descriptive statistics of the real data. Since the real insurance portfolio is confidential, we also design an algorithm to generate synthetic insurance portfolios of a similar structure as the real one, see Appendix B. Finally, in Appendix Cwe provide sensitivity plots of selected neural networks. Risks 2018,6, 29; doi:10.3390/risks6020029 www.mdpi.com/journal/risks Risks 2018,6, 29 2 of 32 1.1. Description of the Simulation Machine The simulation machine is programmed in the language R. The corresponding .zip -folder can be downloaded from the website: https://people.math.ethz.ch/~wmario/simulation.html This .zip -folder contains all parameters, a file readme.pdf which describes the use of our R-functions, as well as the two R-files Functions.V1 and Simulation.Machine.V1 . The first R-file Functions.V1 contains the two R-functions Feature.Generation and Simulation.Machine . The former is used to generate synthetic insurance portfolios (this is described in more detail in Appendix B) and the latter to simulate the corresponding individual claims histories (this is described in the main body of this manuscript). The R-file Simulation.Machine.V1 demonstrates the use of these two R-functions, also providing a short chain-ladder claims reserving analysis. 1.2. Procedure of Developing the Simulation Machine In recent years, neural networks have become increasingly popular in all fields of machine learning. They have proved to be very powerful tools in classification and regression problems. Their drawbacks are that they are rather difficult to calibrate and, once calibrated, they act almost like black boxes between inputs and outputs. Of course, this is a major disadvantage in interpretation and getting deeper insight. However, the missing interpretation is not necessarily a disadvantage in our project because it implies—in back-testing other methods—that the true data generating mechanism cannot easily be guessed. To construct our individual claims history simulation machine, we design a neural network architecture. This architecture is calibrated to real insurance data consisting of n= 9,977,298 individual claims that have occurred between 1994 and 2005. For each of these individual claims, we have full information of 12 years of claims development as well as the relevant feature information. Together with a portfolio generating algorithm (see Appendix B), one can then use the calibrated simulation machine to simulate as many individual claims development histories as desired. 1.3. The Chosen Data The chosen data has been preprocessed correcting for wrong entries—for instance, an accident date that is bigger than the reporting date, etc. Moreover, we have dropped claims with missing feature components—for instance, if the age of the injured was missing. However, this was a negligible number of claims that we had to drop, and this does not distort the general calibration. The final (cleaned) data set consists of n= 9,977,298 individual claims histories. The following feature information is available for each individual claim: •the claims number ClNr, which serves as a distinct claims identifier; •the line of business LoB, which is categorical with labels in {1, . . . , 4}; • the claims code cc , which is categorical with labels in { 1, . . . , 53 } and denotes the labor sector of the injured; •the accident year AY, which is in {1994, . . . , 2005}; •the accident quarter AQ, which is in {1, . . . , 4}; •the age of the injured age (in 5 years age buckets), which is in {15, 20, . . . , 70}; • the injured part inj_part , which is categorical with labels in { 10, . . . , 99 } and denotes the part of the body injured; •the reporting year RY, which is in {1994, . . . , 2016}. Not all values in { 10, . . . , 99 } are needed for the labeling of the categorical classes of the feature component inj_part . In fact, only 46 different values are attained, but for simplicity, we have decided to keep the original labeling received from the insurance company. 46 different labels may still seem to Risks 2018,6, 29 3 of 32 be a lot and a preliminary classification could allow to reduce this number, here we refrain from doing so because each label has sufficient volume. For all claims i= 1, . . . , n , we are given the individual claims cash flow (C(j) i)0≤j≤11 , where C(j) i is the payment for claim i in calendar year AYi+j —and where AYi denotes the accident year of claim i. Note that we only consider yearly payments, i.e., multiple payments and recovery payments within calendar year AYi+j are aggregated into a single, annual payment C(j) i . This single, annual payment can either be positive or negative, depending on having either more claim payments or more recovery payments in that year. The sum over all yearly payments ∑jC(j) iof a given claim ihas to be non-negative because recoveries cannot exceed payments (this is always the case in the considered data). Remark that our simulation machine will allow for recoveries. Finally, for claims i= 1, . . . , n , we are given the claim status process (I(j) i)0≤j≤11 determining whether claim i is open or closed at the end of each accounting year. More precisely, if I(j) i= 1, claim i is open at the end of accounting year AYi+j , and if I(j) i= 0, claim i is closed at the end of that accounting year. Our simulation machine also allows for re-opening of claims, which is quite common in our real data. More description of the data is given in Appendix A. 2. Design of the Simulation Machine Using Neural Networks In this section we describe the architecture of our individual claims history simulation machine. It consists of eight modeling steps: (1) reporting delay T simulation; (2) payment indicator Z simulation; (3) number of payments K simulation; (4) total claim size Y simulation; (5) number of recovery payments K− simulation; (6) recovery size Y− simulation; (7) cash flow (C(j) i)0≤j≤11 simulation and (8) claim status (I(j) i)0≤j≤11 simulation. Each of these eight modeling steps is based on one or several feed-forward neural networks. We introduce the precise setup of such a neural network in Section 2.1 for the simulation of the reporting delay T . Before, we present a global overview of the architecture of our simulation machine. Afterwards, in Sections 2.1–2.8, each single step is described in detail. To start with, we define the initial feature space X1 consisting of the original six feature components as X1={(LoB,cc,AY,AQ,age,inj_part)}. (1) Observe that we drop the claims number ClNr because it does not have explanatory power. Apart from these six feature values, the only other model-dependent input parameters of our simulation machine are the standard deviations for the total individual claim sizes and the total individual recoveries, see Sections 2.4 and 2.6 below. During the simulation procedure, not all of the subsequent steps (1)–(8) may be necessary—e.g., if we do not have any payments, then there is no need to simulate the claim size or the cash flow pattern. We briefly describe the eight modeling steps (1)–(8). (1) In the first step, we use the initial feature space X1 to model the reporting delay T indicating the annualized difference between the reporting year and the accident year. (2) For the second step, we extend the initial feature space X1 by including the additional information of the reporting delay T, i.e., we set X2={(LoB,cc,AY,AQ,age,inj_part,T)}. (2) We use X2to model the payment indicator Zdetermining whether we have a payment or not. (3) For the third step, we set X3=X2and model the number of (yearly) payments K. (4) In the fourth step, we extend the feature space X3 by including the additional information of the number of payments K, i.e., we set X4={(LoB,cc,AY,AQ,age,inj_part,T,K)}, (3) which is used to model the total individual claim size Y. Risks 2018,6, 29 4 of 32 (5) In the fifth step, we model the number of recovery payments K− . We therefore work on the extended feature space X5={(LoB,cc,AY,AQ,age,inj_part,T,K,Y)}. (4) (6) In the sixth step, we model the total individual recovery Y− . To this end, we set X6=X5 . We understand the total individual claim size Y to be net of recovery Y− . Thus, the total payment from the insurance company to the insured is Y+Y− , paid in K−K− yearly payments. The total recovery from the insured to the insurance company is Y−, paid in K−yearly payments. (7) In the seventh step, the task is to generate the cash flows (C(j) i)0≤j≤11 . Therefore, we have to split the total gross claim amount Y+Y− into K−K− positive payments and the total recovery Y− into K− negative payments and distribute these K payments among the 12 development years. For this modeling step, we use different feature spaces X7a, . . . , X7g, all being a subset of X7=LoB,cc,AY,AQ,age,inj_part,T,K,Y,K−,Y−, (5) see Section 2.7 below for more details. (8) In the last step, we model the claim status process (I(j) i)0≤j≤11 , where we use the feature space X8=LoB,AQ,T,C(j)0≤j≤11. Each of these eight modeling steps (1)–(8) consists of one or even multiple feature-response problems, for which we design neural networks. In the end, the full individual claims history simulation machine consists of 35 neural networks. We are going to describe this neural network architecture in more detail next. We remark that some of these networks are rather similar. Therefore, we present the first neural network in full detail, and for the remaining neural networks we focus on the differences to the previous ones. 2.1. Reporting Delay Modeling To model the reporting delay, we work with the initial feature space X1 given in (1) . Let n1=n= 9,977,298 be the number of individual claims in our data. We consider the (annualized) reporting delays Ti, for i=1, . . . , n1, given by Ti=RYi−AYi∈ T ={0, . . . , 11}, where AYi is the accident year and RYi the reporting year of claim i . For confidentiality reasons, we have only received data on a yearly time scale (with the additional information of the accident quarter AQ). A more accurate modeling would use a finer time scale. The three feature components LoB , cc and inj_part are categorical. For neural network modeling, we need to transform these categorical feature components to continuous ones. This could be done by dummy coding, but we prefer the following version because it leads to less parameters. We replace, for instance, the claims code cc by the sample mean of the reporting delay restricted to the corresponding feature label, i.e., for claims code cc =a, we set a7→ a∗=a∗(a) = ∑n1 i=1Ti1{cci=a} ∑n1 i=11{cci=a}∈R, (6) where cc1 , . . . , ccn1 are the observed claims codes. By slight abuse of notation, we obtain a d= 6 dimensional feature space X1 where we may assume that all feature components of X1 are continuous. Such feature pre-processing as in (6) will be necessary throughout this section for the components LoB , Risks 2018,6, 29 5 of 32 cc and inj_part : we just replace Ti in (6) by the respective response variable. Note that from now on this will be done without any further reference. The above procedure equips us with the data D1={(x1,T1), . . . , (xn1,Tn1)}, with x1 , . . . , xn1∈ X1 being the observed features and T1 , . . . , Tn1∈ T the observed responses. For an insurance claim with feature x∈ X1 , the corresponding reporting delay T(x) is modeled by a categorical distribution P[T(x) = t]=πt(x), for t∈ T . This requires that we model probability functions of the form πt:X1→[0, 1],x7→ πt(x), satisfying normalization ∑t∈T πt(x) = 1, for all x∈ X1 . We design a neural network for the modeling of these probability functions and we estimate the corresponding network parameters from the observations D1. We choose a classical feed-forward neural network with multiple layers. Each layer consists of several neurons, and weights connect all neurons of a given layer to all neurons of the next layer. Moreover, we use a non-linear activation function to pass the signals from one layer to the next. The first layer—consisting of the components x1 , . . . , xd of a feature x= (x1 , . . . , xd)∈ X1 —is called input layer (blue circles in Figure 1). In our case, we have d= 6 neurons in this input layer. The last layer is called output layer (red circles in Figure 1) and it contains the categorical probabilities π0(x) , . . . , π11(x) . In between these two layers, we choose two hidden layers having q1 and q2 hidden neurons, respectively (black circles in Figure 1with q1=11 and q2=15). input layer hidden layer hidden layer output layer Figure 1. Deep neural network with two hidden layers: the first column (blue circles) illustrates the d= 6 dimensional feature vector x (input layer), the second column gives the first hidden layer with q1= 11 neurons, the third column gives the second hidden layer with q2= 15 neurons and the fourth column gives the output layer (red circle) with 12 neurons. More formally, we choose the q1hidden neurons z(1) 1, . . . , z(1) q1in the first hidden layer as follows z(1) j=z(1) j(x) = φ w(1) j,0 + d ∑ l=1 w(1) j,lxl!, for all j=1, . . . , q1, Risks 2018,6, 29 6 of 32 for given weights (w(m) j,l)j,l,mand for the hyperbolic tangent activation function φ(x) = tanh(x). This is a centered version of the sigmoid activation function, with range (− 1, 1 ) . Moreover, we have φ0=1−φ2, which is a useful property in the gradient descent method described in Section 3, below. The activation is then propagated in an analogous fashion to the q2 hidden neurons z(2) 1 , . . . , z(2) q2 in the second hidden layer, that is, we set z(2) j=z(2) j(x) = φ w(2) j,0 + q1 ∑ l=1 w(2) j,lz(1) l(x)!, for all j=1, . . . , q2. For the 12 neurons π0(x) , . . . , π11(x) in the output layer, we use the multinomial logistic regression assumption πt(x) = exp {µt(x)} ∑s∈T exp {µs(x)}, for all t∈ T , (7) with regression functions x7→ µt(x)for all t∈ T given by µt(x) = β(t) 0+ q2 ∑ j=1 β(t) jz(2) j(x), for given weights (β(t) j)j,t. We define the network parameter αof all involved parameters by α=w(1) 1,0 , . . . , w(2) q2,q1,β(0) 0, . . . , β(11) q20∈Rq1(d+1)+q2(q1+1)+12(q2+1). The classification model for the tuples (x , T(x))x∈X1 is now fully defined and there remains the calibration of the network parameter α and the choice of the hyperparameters q1 and q2 . Assume for the moment that q1 and q2 are given. In order to fit α to our data D1 , we aim to minimize a given loss function α7→ L(α) . Therefore, we assume that (x1 , T1) , . . . , (xn1 , Tn1) are drawn independently from the joint distribution of (x , T(x)) . The corresponding deviance statistics loss function of the categorical distribution of our data D1is then given by L(α) = LD1(α) = −2 log n1 ∏ i=1 ∑ t∈T 1{Ti=t}πt(xi)!=−2 n1 ∑ i=1 ∑ t∈T 1{Ti=t}log πt(xi). (8) The optimal network parameter α is found by minimizing this deviance statistics loss function. We come back to this problem in Section 3.2.1, below. Since for different hyperparameters q1 and q2 we get different network structures, every pair (q1 , q2) corresponds to a separate model. The choice of appropriate hyperparameters q1and q2is discussed in Section 3.3, below. After the calibration of q1 , q2 and α to our data D1 , we can simulate the reporting delay T(x) of a claim with given feature x∈ X1 by using the resulting categorical distribution given by (7) . This simulated value will then allow us to go to the next modeling step (2), see (2). We close this first part with the following remark: Our choice to work with two hidden layers may seem arbitrary since we could also have chosen more hidden layers or just one of them. From a theoretical point of view, one hidden layer would be sufficient to approximate a vast collection of regression functions to any desired degree of accuracy, provided that we have sufficiently many hidden neurons in that layer, see Cybenko (1989) and Hornik et al. (1989). However, these models with large-scale numbers of hidden neurons are known to be difficult to calibrate, and it is often Risks 2018,6, 29 7 of 32 more efficient to use fewer neurons but more hidden layers to get an appropriate complexity in the regression function. 2.2. Payment Indicator Modeling In our real data, we observe that roughly 29% of all claims can be settled without any payment. For this reason, we model the claim sizes by compound distributions. First, we model a payment indicator Z that determines whether we have a payment or not. Then, conditionally on having a payment, we determine the exact number of payments K . Finally, we model the total individual claim size Yfor claims with at least one payment. In order to model the payment indicator, we work with the d= 7 dimensional feature space X2 introduced in (2) . Let n2=n1 and x1 , . . . , xn2∈ X2 be the observed features, where this time the reporting delay T is also included. For all i= 1, . . . , n2 , we define the number of payments Ki and the payments indicator Ziby Ki= 11 ∑ j=0 1nC(j) i6=0oand Zi=1{Ki>0}. (9) This provides us with the data D2={(x1,Z1), . . . , (xn2,Zn2)}. For a claim with feature x= (x1 , . . . , xd)∈ X2 , the corresponding payment indicator Z(x) is a Bernoulli random variable with P[Z(x) = 1]=π(x), for a given (but unknown) probability function π:X2→[0, 1],x7→ π(x). Note that this Bernoulli model is a special case of the categorical model of Section 2.1. Therefore, it can be calibrated completely analogously, as described above. However, we emphasize that instead of working with two probability functions π0 and π1 for the two categories { 0, 1 } , we set π(·) = π1(·) , which implies 1 −π(·) = 1 −π1(·) = π0(·) . Moreover, the multinomial probabilities (7) simplify to the binomial case π(x) = exp{µ1(x)} exp{µ0(x)}+exp{µ1(x)}=1 1+exp{−(µ1(x)−µ0(x))}=1 1+exp{−µ(x)}, with regression function µ:X2→R,x7→ µ(x) = β0+ q2 ∑ j=1 βjz(2) j(x)(10) for a neural network with two hidden layers and network parameter αgiven by α=w(1) 1,0 , . . . , w(2) q2,q1,β0, . . . , βq20∈Rq1(d+1)+q2(q1+1)+(q2+1). Finally, the corresponding deviance statistics loss function to be minimized is given by L(α) = LD2(α) = −2 n2 ∑ i=1 Zilog π(xi) + (1−Zi)log(1−π(xi)). (11) Risks 2018,6, 29 8 of 32 From this calibrated model, we simulate the payment indicator Z(x) , which then allows us to go to the next modeling step. If this indicator is equal to one, we move to step (3), see Section 2.3; if this indicator is equal to zero, we directly go to step (8), see Section 2.8. 2.3. Number of Payments Modeling We use the d= 7 dimensional feature space X3=X2 to model the number of payments, conditioned on the event that the payment indicator Z is equal to one. We define n3≤n2 to be the number of claims with payment indicator equal to one and order the claims appropriately in i such that Zi= 1 for all i= 1, . . . , n3 . Then, we define the number of payments Ki as in (9) , for all i= 1, . . . , n3 . This gives us the data D3={(x1,K1), . . . , (xn3,Kn3)}. For a claim with feature x∈ X3 and payment indicator Z= 1, we could now proceed as in Section 2.1 in order to model the number of payments K(x) . However, the claims with Ki= 1 are so dominant in the data that a good calibration of the categorical model (7) becomes difficult. For this reason, we choose a different approach: in a first step, we model the events {K(x) = 1 } and {K(x)> 1 } , conditioned on {Z= 1 } , and, in a second step, we consider the conditional distribution of K(x) , given K(x)> 1. In particular, in the first step we have a Bernoulli classification problem that is modeled completely analogously to Section 2.2, only replacing the data D2by D3a=nx1,1{K1=1}, . . . , xn3,1{Kn3=1}o. The case K(x)> 1 is then modeled analogously to the categorical case of Section 2.1, with 11 categories and data D3b⊂ D3only considering the claims with more than one payment. The simulation of the number of payments K(x) for a claim with feature x∈ X3 , reporting delay T and payment indicator Z= 1 needs more care than the corresponding task in Section 2.1: here we have the restriction T+K(x)≤ 12. If T= 11, then we automatically need to have K(x) = 1. For T< 11 and if the first neural network leads to 1{K(x)=1}= 0, then the categorical conditional distribution for K(x) , given K(x)> 1, can only take the values k∈ { 2, . . . , 12 −T} . For this reason, instead of using the original conditional probabilities π2(x) , . . . , π12(x) resulting from the second neural network, we use in that case the modified conditional probabilities π∗ k(x), for k∈ {2, . . . , 12 −T}, given by π∗ k(x) = πk(x) ∑12−T l=2πl(x). (12) 2.4. Total Individual Claim Size Modeling For the modeling of the total individual claim size, we add the number of payments K to the previous feature space and work with X4 given in (3) . Let n4=n3 and consider the same ordering of the claims as in Section 2.3. Then, we define the total individual claim size Yiof claim ias Yi= 11 ∑ j=0 C(j) i>0, for all i= 1, . . . , n4 . In particular, the total individual claim size Yi is always to be understood net of recoveries. This leads us to the data D4={(x1,Y1), . . . , (xn4,Yn4)}. For a claim with feature x∈ X4 and payment indicator Z= 1, we model the total individual claim size Y(x)with a log-normal distribution. We therefore choose a regression function µ:X4→R Risks 2018,6, 29 15 of 32 (c) Proportions of the Recovery Payments if K−=2 In the case of K−= 2 recovery payments, we need to model the proportion of the total individual recovery Y− that is paid in the first recovery payment. For this, we work with the d= 10 dimensional feature space X7g=LoB,cc,AY,AQ,age,inj_part,T,K,Y,Y−. We denote by n7g≤n6 the number of claims with exactly two recovery payments and order the claims appropriately in i such that K− i= 2 for all i= 1, . . . , n7g . Recall that we set K− i= 2 for all claims i with two or more recovery payments, see (15) . Moreover, we add all the amounts of the recovery payments done after the second recovery payment to the second one. Let R− i=min 0≤j≤11 C(j) i<0 denote the time point of the first recovery payment, for all i∈ { 1, . . . , n7g} . Then, the proportion P− i of the total individual recovery Y− ithat is paid in the first recovery payment is given by P− i=−C(R− i) i Y− i , for all i=1, . . . , n7g. This provides us with the data D7g=nx1,P− 1), . . . , (xn7g,P− n7go. The remaining modeling part is then done completely analogously to the second part of the two positive payments case (a) in Section 2.7.2. (d) Cash Flow Modeling Finally, using the three neural network models outlined above, we can simulate the cash flow for a claim with more than two payments and with feature x∈ X7 , see (5) . We illustrate the case K= 6. Note that we only allow for cash flow patterns in e A that are compatible with the reporting delay T . We start by describing the case T= 0. In this case, there is no difficulty and we directly simulate the cash flow pattern R(x)∈e A. This provides us six payments in the time points R(6)(x) = max 0≤j≤11 2j+1≤R(x)and R(k)(x) = max (0≤j<R(k+1)(x)2j+1≤R(x)− 6 ∑ l=k+1 2R(l)(x)+1), for k=1, . . . , 5. For reporting delay T= 1, the set e A of potential cash flow patterns becomes smaller because some of them have to be dropped to remain compatible with T= 1. For this reason, we simulate with probability 1 2 a pattern from e A , and with probability 1 2 the six time points R(1)(x) , . . . , R(6)(x) are drawn in a uniform manner from the remaining possible time points in {T= 1, . . . , 11 } . For T> 1, the potential subset of patterns in e A becomes (almost) empty. For this reason, we simply simulate uniformly from the compatible configurations in {T, . . . , 11}. Having the six time points for the payments, we distinguish the three different cases K−∈ { 0, 1, 2 } : Case K−= 0: we calculate the proportions P(1)(x) , . . . , P(6)(x) according to point (b) above and we receive the cash flow (C(j)(x))0≤j≤11 with C(j)(x) = (P(l)(x)Y, if j=R(l)(x)for some 1 ≤l≤6, 0, else. Risks 2018,6, 29 16 of 32 Case K−= 1: we have five positive payments with proportions P(1)(x) , . . . , P(5)(x) modeled according to point (b) above. This provides the cash flow (C(j)(x))0≤j≤11 with C(j)(x) =      P(l)(x)(Y+Y−), if j=R(l)(x)for some 1 ≤l≤5, −Y−, if j=R(6)(x), 0, else. Case K−= 2:we have four positive payments with proportions P(1)(x) , . . . , P(4)(x) according to point (b) above and two negative payments with proportions P−(x) and 1 −P−(x) according to point (c) above. The time point of the first recovery R−(x) is simulated uniformly from the set of time points {R(2)(x) , . . . , R(5)(x)} . Note that the time point R(1)(x) is reserved for the first positive payment and the time point R(6)(x) for the second recovery payment. We write e R(1)(x) , . . . , e R(4)(x) for the time points of the four positive payments. Summarizing, we get the cash flow (C(j)(x))0≤j≤11 with C(j)(x) =          P(l)(x)(Y+Y−), if j=e R(l)(x)for some 1 ≤l≤4, −P−(x)Y−, if j=R−(x), −(1−P−(x)) Y−, if j=R(6)(x), 0, else. Of course, if K= 3 and K−= 2, we do not need to simulate the proportions of the positive payments, as there is only one positive payment, which occurs in the beginning. Similarly, if K= 12, we do not need to simulate the time points of the payments, since there is a payment in every development year. 2.8. Claim Status Modeling Finally, we design the model for the claim status process which indicates whether a claim is open or closed at the end of each accounting year. This process modeling will also allow for re-opening. Similarly to the payments, we do not model the status of a claim or its changes within an accounting year, but only focus on its status at the end of each accounting year. The modeling procedure of the claim status uses two neural networks, which are described below. We remark that the closing date information was of lower quality in our data set compared to all other information. For instance, some of the dates have been modified retrospectively which, of course, destroys the time series aspect. For this reason, we have decided to model this process in a more crude form, however, still capturing predictive power. 2.8.1. Re-Opening Indicator We start by modeling the occurrence of a re-opening, i.e., whether a claim gets re-opened after having been closed at an earlier date. We use the d=15 dimensional feature space X8a=LoB,AQ,T,e C(j)0≤j≤11, (21) where we do not consider the exact payment amounts, but the simplified version e C(j)=     −1 2, if C(j)=0, 0, if C(j)6=0 and C(j)≤1, 000, 1 2, if C(j)>1, 000, (22) for all j= 0, . . . , 11. Let n8a≤n denote the number of claims i for which we have the full information (I(j) i)Ti≤j≤11 . For the ease of data processing, we set I(j) i= 1 for all development years before claims reporting Ti. Then, we can define the re-opening indicator Vias Risks 2018,6, 29 17 of 32 Vi=   1, if ∑11 j=11nI(j) i−I(j−1) i=1o≥1, 0, else, for all i= 1, . . . , n8a . In particular, if Vi= 1, then claim i has at least one re-opening, and if Vi= 0, then claim ihas not been re-opened. This leads us to the data D8a={(x1,V1), . . . , (xn8a,Vn8a)}, where x1, . . . , xn8a∈ X8a. For a given feature x∈ X8a, the corresponding re-opening indicator V(x)is a Bernoulli random variable. Thus, model calibration is done analogously to Section 2.2 with, however, a neural network with only one hidden layer. 2.8.2. Closing Delay Indicator for Claims without a Re-Opening For claims without a re-opening, we model the closing delay indicator determining whether the closing occurs in the same year as the last payment or if the closing occurs later. In case of no payments ( Zi= 0), we replace the year of the last payment by the reporting year. We use the same d= 15 dimensional feature space as for the re-opening indicator and set X8b=X8a , see (21) . Let n8b≤n8a be the number of claims without a re-opening and order them appropriately in i such that Vi= 0 for all i=1, . . . , n8b. Then, we define the closing delay indicator Wias Wi=       1, if Zi=1 and max n0≤j≤11 I(j) i=1o≥max n0≤j≤11 C(j) i6=0o, 1, if Zi=0 and max n0≤j≤11 I(j) i=1o≥Ti, 0, else, for all i= 1, . . . , n8b . Hence, we have Wi= 1 if the closing occurs in a later year compared to the year of the last payment (or in a later year compared to the claims reporting year in case there is no payment) and Wi=0 otherwise. This leads us to the data D8b={(x1,W1), . . . , (xn8b,Wn8b)}. For a claim with feature x∈ X8b , the corresponding closing delay indicator W(x) is again a Bernoulli random variable. Therefore, model calibration is done analogously to Section 2.2. Similarly as for the re-opening indicator, we use a neural network with only one hidden layer. 2.8.3. Simulation of the Claim Status Based on the feature space X8a , we first simulate the re-opening indicator V(x) leading to the two cases (i) and (ii) described below. Note that before a claim is reported—for ease of data processing—we simply set its status to open (this has no further relevance). (i) Case V(x) = 0 (without re-opening): for the given feature x∈ X8a , we calculate the closing delay probability π(x) = P[W(x) = 1 ] using the neural network of Section 2.8.2. The closing delay B(x)is then sampled from a categorical distribution on {0, . . . , 12}with probabilities P[B(x) = 0] = 1−π(x), P[B(x) = 1] = 9 10 π(x), P[B(x) = k] = 1 10 1 11 π(x), for k=2, . . . , 12. The resulting closing delay B(x)∈ {0, . . . , 12}is added to the year of the last payment (or to the reporting year if there is no payment). If this sum exceeds the value 11, the claim is still open at the end of the last modeled development year. This provides the claim status process (I(j)(x))0≤j≤11 with Risks 2018,6, 29 18 of 32 I(j)(x) =      1, if Z=1 and j<B(x) + max n0≤j≤11 C(j)6=0o, 1, if Z=0 and j<B(x) + T, 0, else. (ii) Case V(x) = 1 (with re-opening): if we have at least one payment for the considered claim, then the first settlement time B1(x)is simulated from a uniform distribution on the set nT, . . . , max n0≤j≤11 C(j)6=0oo. The second settlement time B2(x)is simulated from a uniform distribution on the set nmax n0≤j≤11 C(j)6=0o+2, . . . , 13o. In particular, the first settlement arrives between the reporting year and the year of the last payment. Then, the claim gets re-opened in the year following the first settlement. The second settlement, if there is one, arrives between two years after the year of the last payment and the last modeled development year. In case the second settlement arrives after the last modeled development year, we simply cannot observe it and the claim is still open at the end of the last modeled development year. In case the first settlement happens in the last modeled development year, we do not even observe the re-opening. If the claim does not have any payment, we set B1(x) = T for the first settlement time. In particular, the claim gets closed for the first time in the same year as it is reported. The second settlement time B2(x)is simulated from a uniform distribution on the set {T+2, . . . , 13}. This leads to the claim status process (I(j)(x))0≤j≤11 with I(j)(x) = (1, if j<B1(x)or B1(x)<j<B2(x), 0, if j=B1(x)or j≥B2(x). 3. Model Calibration Using Momentum-Based Gradient Descent In Section 2we have introduced several neural networks that need to be calibrated to the data. This calibration involves the choice of the numbers of hidden neurons q1 and q2 as well as the choice of the corresponding network parameter α . We first focus on the network parameter α for given q1 and q2. 3.1. Gradient Descent Methods State-of-the-art for finding the optimal network parameter α w.r.t. a given differentiable loss function α7→ L(α) is the gradient descent method (GDM). The GDM locally improves the loss in an iterative way. Consider the Taylor approximation of Laround α, then L(eα) = L(α) + (∇αL(α))0(eα−α)+o(keα−αk), as keα−αk → 0. The locally optimal move points into the direction of the negative gradient −∇αL(α) . If we choose a learning rate $>0 into that direction, we obtain a local loss decrease L(α−$∇αL(α))≈ L(α)−$k∇αL(α)k2, (23) for $small. Iterative application of these locally optimal moves—with tempered learning rates—will converge ideally to the (local) minimum of the loss function. Note that (a) it is possible to end up in saddle points; (b) different starting points of this algorithm should be explored to see whether we converge to different (local) minima resp. saddle points and (c) the speed of convergence should be Risks 2018,6, 29 19 of 32 fine-tuned. An improved version of the GDM is the so-called momentum-based GDM introduced in Rumelhart et al. (1986). Consider a velocity vector v with the same dimensions as α and initialize v=0 , corresponding to zero velocity in the beginning. Then, in every iteration step of the GDM, we are building up velocity to achieve a faster convergence. In formulas, this provides v←µv−$∇αL(α), α←α+v, where µ∈[ 0, 1 ] is the momentum coefficient controlling how fast velocity is built up. By choosing µ= 0, we get the original GDM without a velocity vector, see (23) . Fine-tuning 0 <µ≤ 1 may lead to faster convergence. We refer to the relevant literature for more on this topic. 3.2. Gradients of the Loss Functions Involved In Section 2we have met three different model types of neural networks: •categorical case with more than two categorical classes; •Bernoulli case with exactly two categorical classes; •log-normal case. In order to apply the momentum-based GDMs, we need to calculate the gradients of the corresponding loss functions of these three model types. As illustrations, we choose the reporting delay T for the categorical case, the payment indicator Z for the Bernoulli case and the total individual claim size Yfor the log-normal case. 3.2.1. Categorical Case (with More than Two Categorical Classes) The loss function α7→ L(α) for the modeling of the reporting delay T is given in (8) . The gradient ∇αL(α)can be calculated as ∇αL(α) = −2 n1 ∑ i=1 ∑ t∈T 1{Ti=t}∇αlog πt(xi) = −2 n1 ∑ i=1 ∑ t∈T 1{Ti=t} 1 πt(xi)∇απt(xi). We have for the last gradients ∇απt(xi) = ∇α exp {µt(xi)} ∑s∈T exp {µs(xi)}=πt(xi) ∇αµt(xi)−∑ s∈T πs(xi)∇αµs(xi)!, for all t∈ T and i=1, . . . , n1. Collecting all terms, we conclude ∇αL(α) = −2 n1 ∑ i=1 ∑ t∈T 1{Ti=t} ∇αµt(xi)−∑ s∈T πs(xi)∇αµs(xi)! =−2 n1 ∑ i=1 ∑ t∈T 1{Ti=t}−πt(xi)∇αµt(xi). There remains to calculate the gradients ∇αµt(xi) , for all t∈ T and i= 1, . . . , n1 . This is done using the back-propagation algorithm, which in today’s form goes back to Werbos (1982). 3.2.2. Bernoulli Case (Two Categorical Classes) We calculate the gradient ∇αL(α) for the modeling of the payment indicator Z with corresponding loss function L(α)given in (11). We get as in the categorical case above ∇αL(α) = −2 n2 ∑ i=1Zi π(xi)−1−Zi 1−π(xi)∇απ(xi), Risks 2018,6, 29 20 of 32 with gradient ∇απ(xi) = exp {−µ(xi)} (1+exp {−µ(xi)})2∇αµ(xi) = π(xi)(1−π(xi))∇αµ(xi), for all i=1, . . . , n2. Collecting all terms, we obtain ∇αL(α) = −2 n2 ∑ i=1 (Zi−π(xi))∇αµ(xi). We again apply back-propagation to calculate the gradient ∇αµ(xi), for all i=1, . . . , n2. 3.2.3. Log-Normal Case Finally, the loss function L(α) for the modeling of the total individual claim size Y is given in (14) . Hence, for the gradient ∇αL(α), we have ∇αL(α) = n4 ∑ i=1∇α(log Yi−µ(xi))2=−2 n4 ∑ i=1 (log Yi−µ(xi))∇αµ(xi), where the last gradient ∇αµ(xi), for all i=1, . . . , n4, is again calculated using back-propagation. 3.3. Choice of the Numbers of Hidden Neurons For each modeling step of our simulation machine, we still need to determine the optimal neural network in terms of the numbers q1 and q2 of hidden neurons. These hyperparameters are determined by splitting the original data set into a training set and a validation set, where for each calibration we choose at random 90% of the data for the training set. The training set is then used to fit the models for the different choices of hyperparameters q1 and q2 by minimizing the corresponding (training) in-sample losses of the functions α7→ L(α) . This is done as described in the previous sections—for given q1 and q2 . The hyperparameter choices q1 and q2 —and model choices, respectively—are then done by choosing the model with the smallest (validation) out-of-sample loss on the validation set. 4. Chain-Ladder Analysis In this section we use the calibrated stochastic simulation machine to perform a small claims reserving analysis. We generate data from the simulation machine and compare it to the real data. For both data sets, we analyze the resulting claims reporting patterns and the corresponding claims cash flow patterns. For claims reportings, we separate the individual claims i= 1, . . . , n by accident year AY ∈ { 1994, . . . , 2005 } and reporting delays T∈ { 0, . . . , 11 } . For claims cash flows, we separate the individual claims i= 1, . . . , n again by accident year AY ∈ { 1994, . . . , 2005 } and aggregate the corresponding payments over the development delays j= 0, . . . , 11. The reported claims and the claims payments that are available by the end of accounting year 2005 then provide the so-called upper claims reserving triangles. These triangles of reported claims of real and simulated data are shown in Tables 2and 3, the triangles of cumulative claims payments of real and simulated data are given in Tables 4and 5. At a first glance, these triangles show that the simulated data looks very similar to the real data, with a slightly bigger similarity for claims reportings than for claims cash flows. These data sets can be used to perform a chain-ladder (CL) claims reserving analysis. We therefore use Mack’s chain-ladder model, for details we refer to Mack (1993). We calculate the chain-ladder reserves for both the real and the simulated data, and we also calculate Mack’s square-rooted conditional mean square error of prediction √msep. We start the analysis on the claims reportings. Using the chain-ladder method, we predict the number of incurred but not yet reported (IBNYR) claims. These are the predicted numbers of late reported claims in the lower triangles in Tables 2and 3. The resulting predictions are provided in the Risks 2018,6, 29 21 of 32 2nd and 5th columns of Table 6. We observe a high similarity between the results on the real and the simulated data. In particular, for all the individual accident years, the chain-ladder predicted numbers of IBNYR claims of the real data and the simulated data are very close to each other. Aggregating over all accident years, the chain-ladder predicted number of the total IBNYR claims is only 0.2% higher for the simulated data compared to the real data. This similarity largely carries over to the prediction uncertainty analysis illustrated by the columns √msep in Table 6. Indeed, comparing the real and the simulated data, we see that √msep is of similar magnitude for most accident years. Only for the accident years 2003 and 2004 it seems notably higher for the real data. From this, we conclude that, at least from a chain-ladder reserving point of view, our stochastic simulation machine provides very reasonable claims reporting patterns. Finally, Table 7shows the results of the chain-ladder analysis for claims payments. Columns 2 and 5 of that table provide the chain-ladder reserves. These are the payment predictions for the cash flows paid after accounting year 2005 and complete the lower triangles in Tables 4and 5. Also here we see high similarities between the real data and the simulated data analysis: the corresponding total chain-ladder reserves as well as the corresponding reserves for most of the individual accident years are rather close to each other. In particular, the total chain-ladder reserves are only 1.2% higher for the simulated data. We only observe slightly shorter cash flow patterns in the simulated data, which partially carries over to the prediction uncertainties illustrated by the columns √msep in Table 7. Table 2. Triangle of reported claims of the real data. Accident Reporting Delay T Year AY 0 1 2 3 4 5 6 7 8 9 10 11 1994 861,899 59,056 1540 460 230 154 84 56 50 28 32 12 1995 850,297 64,733 1568 562 216 124 94 62 44 34 32 1996 781,875 61,465 1742 414 252 153 76 62 38 22 1997 756,147 59,269 1466 496 210 147 54 48 40 1998 753,552 60,249 1660 530 248 136 98 44 1999 754,992 59,690 1625 468 208 100 68 2000 766,684 61,120 1274 320 136 88 2001 758,443 61,449 1024 286 90 2002 745,125 55,246 876 200 2003 757,843 53,272 956 2004 733,785 51,742 2005 730,978 Table 3. Triangle of reported claims of the simulated data. Accident Reporting Delay T Year AY 0 1 2 3 4 5 6 7 8 9 10 11 1994 860,337 60,143 1837 553 256 155 101 70 57 35 36 21 1995 851,877 62,924 1776 499 263 164 80 58 52 29 36 1996 783,006 60,511 1557 477 182 122 87 49 53 39 1997 756,015 59,556 1434 399 169 129 77 42 42 1998 752,454 61,913 1422 380 164 100 65 37 1999 753,635 61,552 1279 362 153 116 49 2000 768,180 59,636 1215 338 150 71 2001 759,501 60,131 1124 304 132 2002 744,478 55,577 1012 270 2003 757,635 53,352 937 2004 732,884 52,586 2005 731,357 Risks 2018,6, 29 22 of 32 Table 4. Triangle of cumulative claims payments (in 10,000 CHF) of the real data. Accident Development Delay j Year AY 0 1 2 3 4 5 6 7 8 9 10 11 1994 78,433 120,396 130,167 134,749 137,143 138,798 139,994 141,052 142,106 142,913 143,652 144,247 1995 79,372 124,532 135,488 140,338 143,145 144,859 146,408 147,624 148,735 149,686 150,578 1996 71,398 113,335 123,336 128,052 130,779 132,692 134,175 135,359 136,381 137,311 1997 68,600 107,716 117,556 122,331 125,304 127,477 129,037 130,171 131,216 1998 68,055 109,906 120,706 126,443 129,727 131,998 133,745 135,176 1999 71,989 114,344 127,311 134,174 138,075 140,614 142,530 2000 72,225 118,418 131,615 138,216 142,024 144,450 2001 74,891 126,244 141,008 147,923 151,917 2002 78,167 129,105 143,731 150,560 2003 82,668 134,010 148,161 2004 80,630 130,390 2005 82,015 Table 5. Triangle of cumulative claims payments (in 10,000 CHF) of the simulated data. Accident Development Delay j Year AY 0 1 2 3 4 5 6 7 8 9 10 11 1994 80,491 117,807 129,673 135,331 138,591 140,745 142,268 143,394 144,344 145,090 145,694 146,076 1995 79,170 116,943 129,313 135,330 138,785 141,079 142,743 144,045 145,115 146,002 146,760 1996 71,675 107,228 119,578 125,407 128,771 131,009 132,566 133,718 134,657 135,472 1997 68,857 104,291 116,406 122,177 125,495 127,662 129,194 130,388 131,270 1998 67,418 103,415 116,194 122,262 125,673 127,877 129,452 130,588 1999 69,308 107,587 121,381 127,930 131,534 133,924 135,502 2000 73,359 113,266 127,878 134,804 138,806 141,431 2001 73,338 115,626 131,197 138,481 142,591 2002 74,887 117,602 133,850 141,539 2003 81,921 127,461 144,444 2004 81,394 128,864 2005 86,837 Table 6. Chain-ladder predicted numbers of incurred but not yet reported (IBNYR) claims and Mack’s √msep for the real and the simulated data. Accident CL Predicted √msep in % CL Predicted √msep in % Year AY IBNYR Claims IBNYR Claims Real Data Simulated Data 1994 0 0 1995 12 0 0.0% 21 0 0.0% 1996 40 0 0.4% 52 0 0.4% 1997 65 5 8.4% 82 7 8.6% 1998 105 7 6.2% 129 9 6.6% 1999 156 10 6.1% 178 14 7.7% 2000 235 19 8.0% 255 21 8.3% 2001 357 32 9.0% 370 35 9.4% 2002 536 65 12.2% 535 53 9.9% 2003 944 135 14.3% 925 95 10.2% 2004 2201 330 15.0% 2170 265 12.2% 2005 57,734 3542 6.1% 57,789 3410 5.9% total 62,385 3565 5.7% 62,506 3425 5.5% Risks 2018,6, 29 23 of 32 Table 7. Chain-ladder reserves for claims payments (in 10,000 CHF) and Mack’s √msep for the real and the simulated data. Accident CL Reserves √msep in % CL Reserves √msep in % Year AY Real Data Simulated Data 1994 0 0 1995 624 117 18.8% 385 109 28.2% 1996 1337 146 10.9% 991 161 16.2% 1997 2112 168 8.0% 1723 179 10.4% 1998 3224 177 5.5% 2636 187 7.1% 1999 4686 259 5.5% 3943 209 5.3% 2000 6476 394 6.1% 5826 230 3.9% 2001 9275 599 6.5% 8846 290 3.4% 2002 13,049 889 6.8% 12,489 421 3.4% 2003 19,973 1421 7.1% 20,817 837 4.0% 2004 32,532 2394 7.4% 36,647 2050 5.6% 2005 82,706 5039 6.1% 84,175 4968 5.9% total 175,994 6275 3.6% 178,076 5732 3.2% 5. Conclusions We have developed a stochastic simulation machine that generates individual claims histories of non-life insurance claims. This simulation machine is based on neural networks which have been calibrated to real non-life insurance data. The inputs of the simulation machine are a portfolio of non-life insurance claims—for which we want to simulate the corresponding individual claims histories—and the two variance parameters σ2 + (for the total individual claim size, see (13) ) and σ2 − (for the total individual recovery, see Section 2.6). Together with a portfolio generating algorithm, see Appendix B, one can use this simulation machine to simulate as many individual claims histories as desired. In a chain-ladder analysis we have seen that the simulation machine leads to reasonable results, at least from a chain-ladder reserving point of view. Therefore, our simulation machine may serve as a stochastic scenario generator for individual claims histories, which provides a common ground for research in this area, we also refer to the study in Wüthrich (2018b). Acknowledgments: Our greatest thanks go to Suva, Peter Blum and Olivier Steiger for providing data, their insights and for their immense support. Author Contributions: Both authors contributed equally to this work. Conflicts of Interest: The authors declare no conflict of interest. Appendix A. Descriptive Statistics of the Chosen Data Set In this appendi we provide descriptive statistics of the data used to calibrate the individual claims history simulation machine. For confidentiality reasons, we can only show aggregate statistics of the claims portfolio, see Figures A1–A4 below. Risks 2018,6, 29 24 of 32 1 2 3 4 line of business line of business (LoB) number of claims (a) 1 8 15 23 31 40 48 claims code claims code (cc) number of claims (b) 1994 1997 2000 2003 accident year accident year (AY) number of claims (c) 1 2 3 4 accident quarter accident quarter (AQ) number of claims (d) 15 25 35 45 55 65 age of injured age of injured (age) number of claims (e) 10 20 30 37 50 60 71 injured part injured part (inj_part) number of claims (f) Figure A1. Portfolio distributions w.r.t. the features ( a ) LoB ; ( b ) cc ; ( c ) AY ; ( d ) AQ ; ( e ) age and ( f ) inj_part. 0 2 4 6 8 10 reporting delay reporting delay (T) logarithmic number of claims (a) 0 2 4 6 8 10 reporting delay reporting delay (T) average claim size (b) 0 2 4 6 8 10 reporting delay reporting delay (T) average number of payments (c) Figure A2. ( a ) Logarithmic number of claims; ( b ) average claim size and ( c ) average number of payments w.r.t. the reporting delay T; the red lines show the averages. Risks 2018,6, 29 31 of 32 0.2 0.4 0.6 0.8 1.0 first recovery proportion claims code (cc) proportion 1 8 15 23 31 40 48 (a) 0.2 0.4 0.6 0.8 1.0 first recovery proportion injured part (inj_part) proportion 10 20 30 37 50 60 71 (b) 0.2 0.4 0.6 0.8 1.0 first recovery proportion total individual recovery (Y−) proportion 1 10 100 10,000 (c) Figure A14. Proportion P− of the total individual recovery Y− paid in the first recovery payment in the case of K−= 2 recovery payments w.r.t. the features ( a ) cc ; ( b ) inj_part and ( c ) total individual recovery Y−. In Figure A15 we see that claims in lines of business one and four have a higher re-opening probability than claims in lines of business two and three. Moreover, the higher the reporting delay T of a claim, the lower the rate of reopening. Finally, the probability of re-opening heavily depends on the cash flow. In order to not overload the plot, we only show sensitivities w.r.t. the payments C(0) , C(1) , C(2) , C(3) , C(8) , C(10) . Recall that for this neural network, the yearly payments are coded with the values −1 2 , 0 and 1 2 , see (22) . Summarizing, one can say that if we have a payment after the first development year, then the probability of re-opening is quite high. 0.0 0.2 0.4 0.6 re−opening indicator line of business (LoB) probability of re−opening 1 2 3 4 (a) 0.0 0.2 0.4 0.6 re−opening indicator reporting delay (T) probability of re−opening 0246810 (b) 111 0.0 0.2 0.4 0.6 re−opening indicator payments (C0,C1,C2,C3,C8,C10) probability of re−opening 2 2 2 3 33 4 44 5 5 5 6 6 6 1 23 45 6 C0 C1 C2 C3 C8 C10 −0.5 0 0.5 (c) Figure A15. Re-opening indicator V w.r.t. the features ( a ) LoB ; ( b ) reporting delay T and ( c ) yearly payments C(0),C(1),C(2),C(3),C(8),C(10). References Antonio, Katrien, and Richard Plat. 2014. Micro-Level Stochastic Loss Reserving for General Insurance. Scandinavian Actuarial Journal 7: 649–69. Cybenko, George. 1989. Approximation by Superpositions of a Sigmoidal Function. Mathematics of Control, Signals, and Systems (MCSS) 2: 303–14. Hiabu Munir, Carolin Margraff, Maria D. Martínez-Miranda, and Jens P. Nielsen. 2016. The Link between Classical Reserving and Granular Reserving through Double Chain-Ladder and its Extensions. British Actuarial Journal 21: 97–116. Hornik, Kurt, Maxwell Stinchcombe, and Halbert White. 1989. Multilayer Feedforward Networks are Universal Approximators. Neural Networks 2: 359–66. Risks 2018,6, 29 32 of 32 Jessen, Anders H., Thomas Mikosch, and Gennady Samorodnitsky. 2011. Prediction of Outstanding Payments in a Poisson Cluster Model. Scandinavian Actuarial Journal 3: 214–37. Mack, Thomas. 1993. Distribution-Free Calculation of the Standard Error of Chain Ladder Reserve Estimates. ASTIN Bulletin 23: 213–25. Martínez-Miranda, Maria D., Jens P. Nielsen, Richard J. Verrall, and Mario V. Wüthrich. 2015. The Link between Classical Reserving and Granular Reserving through Double Chain-Ladder and its Extensions. Scandinavian Actuarial Journal 5: 383–405. Pigeon, Mathieu, Katrien Antonio, and Michel Denuit. 2013. Individual Loss Reserving with the Multivariate Skew Normal Framework. ASTIN Bulletin 43: 399–428. Rumelhart, David E., Geoffrey E. Hinton, and Ronald J. Williams. 1986. Learning Representations by BackPropagating Errors. Nature 323: 533–36. Taylor, Greg, Gráinne McGuire, and James Sullivan. 2008. Individual Claim Loss Reserving Conditioned by Case Estimates. Annals of Actuarial Science 3: 215–56. Verrall, Richard J., and Mario V. Wüthrich. 2016. Understanding Reporting Delay in General Insurance. Risks 4: 25. Werbos, Paul J. 1982. Applications of Advances in Nonlinear Sensitivity Analysis. In System Modeling and Optimization. Paper Presented at the 10th IFIP Conference, New York City, NY, USA, 31 August–4 September 1981. Edited by Rudolf F. Drenick and Frank Kozin. Berlin and Heidelberg: Springer, pp. 762–70. Wüthrich, Mario V. 2018. Machine Learning in Individual Claims Reserving. To appear in Scandinavian Actuarial Journal 25: 1–16. Wüthrich, Mario V. 2018. Neural Networks Applied to Chain-Ladder Reserving. SSRN Manuscript. doi:10.2139/ssrn.2966126. c 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).