Nagging predictors

Richman, Ronald; Wüthrich, Mario V. Article Nagging predictors Risks Provided in Cooperation with: MDPI – Multidisciplinary Digital Publishing Institute, Basel Suggested Citation: Richman, Ronald; Wüthrich, Mario V. (2020) : Nagging predictors, Risks, ISSN 2227-9091, MDPI, Basel, Vol. 8, Iss. 3, pp. 1-26, https://doi.org/10.3390/risks8030083 This Version is available at: https://hdl.handle.net/10419/258036 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 Nagging Predictors Ronald Richman 1,† and Mario V. Wüthrich 2,*,† 1QED Actuaries and Consultants, 38 Wierda Road West, Sandton 2196, South Africa; [email protected] 2RiskLab, Department of Mathematics, ETH Zurich, 8092 Zurich, Switzerland *Correspondence: [email protected] † These authors contributed equally to this work. Received: 24 June 2020; Accepted: 24 July 2020; Published: 4 August 2020   Abstract: We define the nagging predictor, which, instead of using bootstrapping to produce a series of i.i.d. predictors, exploits the randomness of neural network calibrations to provide a more stable and accurate predictor than is available from a single neural network run. Convergence results for the family of Tweedie’s compound Poisson models, which are usually used for general insurance pricing, are provided. In the context of a French motor third-party liability insurance example, the nagging predictor achieves stability at portfolio level after about 20 runs. At an insurance policy level, we show that for some policies up to 400 neural network runs are required to achieve stability. Since working with 400 neural networks is impractical, we calibrate two meta models to the nagging predictor, one unweighted, and one using the coefficient of variation of the nagging predictor as a weight, finding that these latter meta networks can approximate the nagging predictor well, only with a small loss of accuracy. Keywords: bagging; bootstrap aggregation; neural networks; network aggregation; insurance pricing; regression modeling 1. Introduction Aggregating is a statistical technique that helps to reduce noise and uncertainty in predictors. Typically, it is applied to an i.i.d. sequence of predictors and it is justified theoretically using the law of large numbers. In many practical applications one is not in the comfortable situation of having an i.i.d. sequence of predictors to which aggregating could be applied. For this reason, Breiman (1996) combined b ootstrapping and agg regat ing , called bagging , bootstrapping being used to generate a sequence of randomized predictors and then aggregating them to receive an averaged bootstrap predictor. The title of this paper has been inspired by Breiman (1996), in other words, it is not meant in the sense of grumbling or the like, but we are going to combine n etworks and agg regat ing to receive the nagging predictor. Thereby, we benefit from the fact that typically neural network regression models have infinitely many equally good predictors which are determined by gradient descent algorithms. These equally good predictors depend on the (random) starting point of the gradient descent algorithm, thus, starting the algorithm in i.i.d. seeds for different runs results in an infinite sequence of (equally good) neural network predictors. Moreover, other common neural network training techniques may even lead to more randomness in neural network results, for example, dropout, which relies on randomly setting parts of the network to zero during training to improve performance at test time, see Srivastava et al. (2014) for more details, or the random selection of data to implement stochastic or mini-batch gradient descent methods. On the one hand, this creates difficulties for using neural network models within an insurance pricing context, since the predictive performance of the models measured at portfolio level will vary with each run and the predictions for individual policies will vary Risks 2020,8, 83; doi:10.3390/risks8030083 www.mdpi.com/journal/risks Risks 2020,8, 83 2 of 26 even more, leading to uncertainty about the prices that should ultimately be charged to individuals. On the other hand, having multiple network predictors puts us in the same situation as Breiman (1996) after having received the bootstrap samples, and we can aggregate them, leading to more stable results and enhanced predictive performance. In this paper, we explore statistical properties of nagging predictors at a portfolio and at a policy level, in the context of insurance pricing. We give a short review of related literature. Aggregation of predictors is also known as ensembling. Dietterich (2000a,2000b) discusses ensembling within various machine learning methods such as decision trees and neural networks. As stated in Dietterich (2000b), successful ensembling crucially relies on the fact that one has to be able to construct multiple suitable predictors. This is usually achieved by injecting randomness either into the data or into the fitting algorithm. In essence, this is what we do in this study, however, our main objective is to reduce randomness (in the results of the fitting algorithms) because, from a practical perspective, we need to have uniqueness in insurance pricing, in particular, we prove rates of convergence at which the random element can be diminished. Zhou (2012) and Zhou et al. (2002) study optimal ensembling of multiple predictors. Their proposal leads to a quadratic optimization problem that can be solved with the method of Lagrange. This is related to our task below of finding an optimal meta model using individual uncertainties which we discuss in Section 5.8, below. Related empirical studies on financial data are given in Di Persio and Honchar (2016), du Jardin (2016) and Wang et al. (2011). These studies conclude that ensembling is very beneficial to improve predictive performance which, indeed, is also one of our main findings. Organization of manuscript. In the next section we introduce the framework of neural network regression models, and we discuss their calibration within Tweedie’s compound Poisson models. In Section 3we give the theoretical foundation of aggregating predictors, in particular, we prove that aggregating leads to more stability in prediction in terms of a central limit theorem. In Section 4we implement aggregation within the framework of neural network predictors providing the nagging predictor. In Section 5we empirically prove the effectiveness of nagging predictors based on a French car insurance data set. Since nagging predictors involve simultaneously dealing with multiple networks, we provide a more practical meta network that approximates the nagging predictor at a minimal loss of accuracy. 2. Feed-Forward Neural Network Regression Models These days neural networks are state-of-the-art for performing complex regression and classification tasks, particularly on unstructured data such as images or text. For a general introduction to neural networks we refer to LeCun et al. (2015), Goodfellow et al. (2016) and the references therein. In the present work, we follow the terminology and notation of Wüthrich (2019). We design feed-forward neural network regression models, referred to in short as networks, to predict insurance claims. In the spirit of Wüthrich (2019), we understand these networks as extensions of generalized linear models (GLMs), the latter being introduced by Nelder and Wedderburn (1972). 2.1. Generic Definition of Feed-Forward Neural Networks Assume we have independent observations (Yi , xi , vi) , i= 1, . . . , n ; the variables Yi describe the responses (here: insurance claims), xi∈ X ⊂ Rq0 describe the real-valued feature information (also known as covariates, explanatory variables, independent variables or predictors), and vi> 0 are known exposures. The main goal is to find a regression functional µ(·) that appropriately describes the expected insurance claims E[Yi]as a function of the feature information xi, i.e., X → R,xi7→ µ(xi) = E[Yi]. Risks 2020,8, 83 3 of 26 Since, typically, the true regression function µ is unknown we approximate it by a network regression function. Under a suitable link function choice g(·) , we assume that E[Yi] can be described by the following network of depth d∈N xi7→ g(E[Yi])=Dβ(d+1),z(d)◦···◦z(1)(xi)E, (1) where h· , ·i is the scalar product in Euclidean space Rqd+1 , the operation ◦ gives a composition of hidden network layers z(m) , 1 ≤m≤d , of dimensions qm+ 1 ∈N , and β(d+1)∈Rqd+1 is the readout parameter. For a given non-linear activation function φ:R→R , the m -th hidden network layer is a mapping z(m):{1}×Rqm−1→ {1}×Rqm,x7→ z(m)(x) = 1, z(m) 1(x), . . . , z(m) qm(x)>, with hidden neurons z(m) j(·), 1 ≤j≤qm, being described by ridge functions x7→ z(m) j(z) = φDβ(m) j,xE, for network parameters β(m) j∈Rqm−1+1 . This network regression function has a network parameter β= (β(1) 1, . . . , β(q) qd,β(d+1))∈Rrof dimension r=∑d+1 m=1(qm−1+1)qm. Based on so-called universality theorems, see, for example, Cybenko (1989) and Hornik et al. (1989) , we know that networks are very flexible and can approximate any continuous and compactly supported regression function arbitrarily well if one allows for a sufficiently complex network architecture. In practical applications this means that one has to choose sufficiently many hidden layers with sufficiently many hidden neurons to be assured of adequate approximation capacity, and then one can work with this network architecture as a replacement of the unknown regression function. This complex network is calibrated to (learning) data using the gradient descent algorithm. To prevent the network from over-fitting to the learning data, one exercises early stopping of this calibration algorithm by applying a given stopping rule. This early stopping of the algorithm before convergence typically implies that one cannot expect to receive a “unique best” model, in fact, the early stopped calibration depends on the starting point of the gradient descent algorithm and, usually, since each run has a different starting point, one receives a different solution each time early-stopped gradient descent is run. As described in Section 3.3.1 of Wüthrich (2019), this exactly results in the problem of having infinitely many equally good models for a fixed stopping rule, and it is not clear which particular solution should be chosen, say, for insurance pricing. This is the main point on which we elaborate in this work, and we come back to this discussion after Equation (7), below. 2.2. Tweedie’s Compound Poisson Model We assume that the response Yi in Equation (1) belongs to the family of Tweedie’s compound Poisson (CP) models, see Tweedie (1984), having a density of the form Yi∼f(y;θi,vi/ϕ,p) = exp yθi−κp(θi) ϕ/vi +ap(y;vi/ϕ), (2) Risks 2020,8, 83 4 of 26 with p∈[1, 2]being the power variance (hyper-)parameter, and vi>0 is a given exposure (weight, volume), ϕ>0 is the dispersion parameter, θi∈Θpis the canonical parameter in the effective domain Θp, κp:Θp→Ris a cumulant function of the form Equation (3), below, ap(·;·)is the normalization, not depending on the canonical parameter θi. Tweedie’s CP models belong to the exponential dispersion family (EDF), the latter being more general because it allows for more general cumulant functions κ , we refer to Jørgensen (1986,1987). We focus here on Tweedie’s CP family and not on the entire EDF because for certain calculations we need explicit properties of cumulant functions. However, similar results can be derived for other members of the EDF. We mention that the effective domain Θp⊂R is a convex set, and that the cumulant function κp is smooth and strictly convex in the interior of the effective domain, having the following explicit form κp(θ) =      exp(θ)for p=1, 1 2−p((1−p)θ)2−p 1−pfor p∈(1, 2), −log(−θ)for p=2. (3) p= 1 is the Poisson model (with ϕ= 1), p= 2 is the gamma model, and for p∈( 1, 2 ) we receive compound Poisson models with i.i.d. gamma claim sizes having shape parameter γ= ( 2 −p)/(p− 1 )∈( 0, ∞) , see Jørgensen and de Souza (1994) and Smyth and Jørgensen (2002). The first two moments are given by µi def. =E[Yi] = κ0p(θi)and Var(Yi) = ϕ vi κ00 p(θi) = ϕ vi Vp(µi), (4) with power variance function Vp(µ) = µpfor p∈[1, 2]. Remark 1. In this paper we work under Tweedie’s CP models characterized by cumulant functions of the form Equation (3) because we need certain (explicit) properties of κp in the proofs of the statements below. For regression modeling one usually starts from a bigger class of models, namely, the EDF, which only requires that the cumulant function is smooth and strictly convex on the interior of the corresponding effective domain Θ , see Jørgensen (1986,1987). A sub-family of the EDF is the so-called Tweedie’s distributions which have a cumulant function that allows for a power mean-variance relationship Equation (4) for any p∈R\( 0, 1 ) , see Table 1 in Jørgensen (1987), for instance, p= 0is the Gaussian model, p= 1is the Poisson model, p= 2is the gamma model and p= 3is the inverse Gaussian model. As mentioned, the models generated by cumulant function Equation (3) cover the interval ( 1, 2 ) and correspond to compound Poisson models with i.i.d. gamma claim size, we also refer to Delong et al. (2020). Relating Equation (3) to network regression function Equation (1) implies that we aim at modeling the canonical parameter θiof policy iby xi7→ θi=θ(xi) = κ0p−1(µ(xi))=κ0p−1g−1Dβ(d+1),z(d)◦···◦z(1)(xi)E, (5) and (κ0p)−1 is the canonical link, in contrast to the general link function g . If we choose the canonical link for g then (κ0p)−1◦g−1 provides the identity function in Equation (5) . Network parameter β∈Rr is estimated with maximum likelihood estimation (MLE). This is achieved either by maximizing the log-likelihood function or by minimizing the corresponding deviance loss function. We use the framework of the deviance loss function here because it is more closely related to the philosophy Risks 2020,8, 83 5 of 26 of minimizing an objective function in gradient descent methods. The average deviance loss for independent random variables Yiis under Tweedie’s CP model assumption given by L(D;β) = 2 n∑n i=1 vi ϕYiκ0p−1(Yi)−κpκ0p−1(Yi)−Yiθi+κp(θi) =2 n∑n i=1 vi ϕYiκ0p−1(Yi)−κpκ0p−1(Yi)−Yiκ0p−1(µi)+κpκ0p−1(µi), (6) for given data D={(Yi , xi , vi) ; i= 1, . . . , n} . The canonical parameter is given by θi=θ(xi , β) if we understand it as a function defined through Equation (5) . Naturally we have for any mean parameter µiof any policy i δ(Yi,µi)def. =2vi ϕYiκ0p−1(Yi)−κpκ0p−1(Yi)−Yiκ0p−1(µi)+κpκ0p−1(µi) ≥0, (7) because these terms subtract twice the log-likelihood of the µi -parametrized model from its saturated counterpart. δ(Yi,µi)is called the unit deviance of Yiw.r.t. mean parameter µi. The gradient descent algorithm now tries to make the objective function L(D ; β) small through optimizing network parameter β . This is done globally, i.e., simultaneously on the entire portfolio i= 1, . . . , n . Naturally, two different parameters β(1)6=β(2) with L(D ; β(1)) = L(D ; β(2)) may provide very different models on an individual policy level i . This is exactly the point raised in the previous section, namely, that early stopping in model calibration w.r.t. to a global objective function L on observations D may provide equally good models on that portfolio level, but they may be very different on an individual policy level (if the gradient descent algorithm has not converged to the same extremal point of the loss function). The goal of this paper is to study such differences on individual policy level. 3. Aggregating Predictors We introduce aggregating in this section. This can be described on one single policy i . Assume that b µi is a predictor for response Yi (and an estimator for mean parameter µi=E[Yi] ). In general, we assume that predictor b µi and response Yi are independent. This independence reflects that we perform an out-of-sample analysis, meaning that the mean parameter has been estimated on a learning data set that is disjoint from Yi , and we aim at performing a generalization analysis by considering the average loss described by the expected unit deviance (subject to existence) E[δ(Yi,b µi)]. (8) We typically assume that the predictor b µi is chosen such that this expected generalization loss is finite. Proposition 1. Choose response Yi∼f(· ; θi , vi/ϕ , p) with power variance parameter p∈[ 1, 2 ] and canonical parameter θi∈Θp . Assume b µi is an unbiased estimator for the mean parameter µi=κ0p(θi) , being independent of Yi , and additionally satisfying e<b µi≤p/(p− 1 )µi , a.s., for some e∈( 0, p/(p− 1 )µi) . We have expected generalization loss E[δ(Yi,b µi)]≥E[δ(Yi,µi)]. We remark that the bounds on b µi ensure that Equation (8) is finite, moreover, the upper bound p/(p− 1 )µi is needed to ensure that the predictors lie in the domain such that the expected unit deviances m7→ E[δ(Yi,m)] are convex functions in m. This is needed in the following proofs. Proof. We calculate the expected generalization loss received by the mean of the unit deviance Risks 2020,8, 83 6 of 26 E[δ(Yi,b µi)]=2vi ϕEYiκ0p−1(Yi)−κpκ0p−1(Yi)−Yiκ0p−1(b µi)+κpκ0p−1(b µi) =E[δ(Yi,µi)] +2vi ϕEYiκ0p−1(µi)−κpκ0p−1(µi)−Yiκ0p−1(b µi)+κpκ0p−1(b µi) =E[δ(Yi,µi)]+2vi ϕhp(µi)−Ehp(b µi), where in the last step we have used independence between Yiand b µiand where we use function m>07→ hp(m) = µiκ0p−1(m)−κpκ0p−1(m) =     µilog(m)−mfor p=1, µim1−p 1−p−m2−p 2−pfor p∈(1, 2), −µi/m−log(m)for p=2. (9) We calculate the second derivative of this function. For p∈[1, 2], it is given by ∂2 ∂m2hp(m) = −pµim−p−1−(1−p)m−p=m1+p[−pµi−(1−p)m]≤0, where for the last inequality we have used that the square bracket is non-negative under our assumptions. This implies that hpis a concave function, and applying Jensen’s inequality we obtain E[δ(Yi,b µi)]=E[δ(Yi,µi)]+2vi ϕhp(µi)−Ehp(b µi) ≥E[δ(Yi,µi)]+2vi ϕhp(µi)−hp(E[b µi])=E[δ(Yi,µi)], where in the last step we have used unbiasedness of estimator b µi. This finishes the proof. Proposition 1tells us that the estimated model b µi has an expected generalization loss Equation (8) which is bounded below by the one of the true model mean µi of Yi . Using aggregating we now try to come as close as possible to this lower bound. Breiman (1996) has analyzed this question in terms of the square loss function, and further results under the square loss function are given in Bühlmann and Yu (2002). We prefer to work with the deviance loss function here because this is the objective function used for fitting in the gradient descent algorithm. For this reason we prove the subsequent results, however, in their deeper nature these results are equivalent to the ones in Breiman (1996) and Bühlmann and Yu (2002). Assume that b µ(j) iare i.i.d. copies of unbiased predictor b µi. We define the aggregated predictor ¯ µ(M) i=1 M M ∑ j=1b µ(j) i. (10) Proposition 2. Assume that b µ(j) i , j≥ 1, are i.i.d. copies of b µi satisfying the assumptions of Proposition 1, and being all independent from Yi. We have for all M ≥1 EhδYi,b µ(1) ii ≥EhδYi,¯ µ(M) ii ≥EhδYi,¯ µ(M+1) ii ≥E[δ(Yi,µi)]. Risks 2020,8, 83 7 of 26 Proof. The last bound is immediately clear because the aggregated predictors themselves fulfill the assumptions of Proposition 1, and henceforth the corresponding statement. Thus, we focus on the inequalities for M≥1. Consider decomposition of the aggregate predictor for M+1 ¯ µ(M+1) i=1 M+1 M+1 ∑ j=1 ¯ µ(−j) i, where ¯ µ(−j) i=1 M M+1 ∑ k=1b µ(k) i1{k6=j}. The predictors ¯ µ(−j) i , j≥ 1, are copies of ¯ µ(M) i , though not independent ones. We have, using function hpdefined in Equation (9), Ehδ(Yi,¯ µ(M) i)i=Ehδ(Yi,¯ µ(M+1) i)i+2vi ϕEhhp¯ µ(M+1) ii−Ehhp¯ µ(M) ii =Ehδ(Yi,¯ µ(M+1) i)i+2vi ϕ E"hp 1 M+1 M+1 ∑ j=1 ¯ µ(−j) i!#−Ehhp¯ µ(M) ii! ≥Ehδ(Yi,¯ µ(M+1) i)i+2vi ϕ E"1 M+1 M+1 ∑ j=1 hp¯ µ(−j) i#−Ehhp¯ µ(M) ii! =Ehδ(Yi,¯ µ(M+1) i)i, where the inequality follows from applying Jensen’s inequality to the concave function hp , and the last identity follows from the fact that ¯ µ(−j) i,j≥1, are copies of ¯ µ(M) i. This finishes the proof. Proposition 2says that aggregation works, i.e., aggregating i.i.d. predictors Equation (10) leads to a monotonically decreasing expected generalization loss. Moreover, notice that the i.i.d. assumption can be relaxed, indeed, it is sufficient that every ¯ µ(−j) i in the above proof has the same distribution as ¯ µ(M) i . This does not require independence between the predictors b µ(j) i , j≥ 1, but exchangeability is sufficient. Proposition 3. Assume that b µ(j) i , j≥ 1, are i.i.d. copies of b µi satisfying the assumptions of Proposition 1, and being all independent from Yi . In the Poisson case p= 1we additionally assume that the sequence of aggregated predictors Equation (8)has a uniform integrable upper bound. We have lim M→∞ EhδYi,¯ µ(M) ii=Elim M→∞δYi,¯ µ(M) i=E[δ(Yi,µi)]. Proof. We have the identity Ehδ(Yi,¯ µ(M) i)i=2vi ϕEYiκ0p−1(Yi)−Eκpκ0p−1(Yi)−Ehhp¯ µ(M) ii, thus, it suffices to consider the last term. The law of large numbers implies a.s. convergence limM→∞hp(¯ µ(M) i) = hp(µi) because we have i.i.d. unbiased predictors b µ(j) i , j≥ 1 (in particular we have consistency). Thus, it suffices to provide a uniform integrable bound and then the claim follows from Lebesgue’s dominated convergence theorem. Note that by assumption we have uniform bounds ¯ µ(M) i∈(e , p/(p− 1 )µi) , a.s., which proves the claim for p∈( 1, 2 ] . Thus, there remains the Poisson case p= 1. In the Poisson case we have h1(m) = µilog(m)−m . The leading term of this function is linear for m→∞ , henceforth, the uniform integrable bound assumption on the sequence Equation (8) provides the proof. Risks 2020,8, 83 8 of 26 The previous statement is based on the law of large numbers. Of course, we can also study a central limit theorem (CLT) that provides asymptotic normality and the rate of convergence. For the aggregated predictors we have convergence in distribution M1/2 ¯ µ(M) i−µi Var(b µi)1/2 =⇒ N(0, 1), as M→∞, (11) noting that, under the Poisson case p= 1 we need to assume, in addition to the assumptions of Proposition 1, that the second moment of b µi is finite. Consider the expected generalization loss function on policy ifor given mean estimate m>0 being independent of Yi m>07→ ¯ δi(m) = E[δ(Yi,m)|m]=2vi ϕEYiκ0p−1(Yi)−Eκpκ0p−1(Yi)−hp(m), where function hp was defined in Equation (9) . Thus, for a CLT of the expected generalization loss it suffices to understand the asymptotic behavior of hp(¯ µ(M) i) , because we have, using the tower property for conditional expectation and using independence between ¯ µ(M) iand Yi, EhδYi,¯ µ(M) ii=EhEhδYi,¯ µ(M) i¯ µ(M) iii=Eh¯ δi¯ µ(M) ii. Proposition 4. Assume that b µ(j) i , j≥ 1, are i.i.d. copies of b µi satisfying the assumptions of Proposition 1, and all being independent from Yi . In the Poisson case p= 1we additionally assume that b µi has finite second moment. We have M1/2 hp(¯ µ(M) i)−hp(µi) h0p(µi)Var(b µi)1/2 =⇒ N(0, 1),as M →∞. This shows how the speed of convergence of aggregated predictors translates to the speed of convergence of deviance loss functions, in particular, from Taylor’s expansion we receive a first order term h0p(µi) . Basically, this is Theorem 5.2 in Lehmann (1983); we provide a proof because it is instructive. Proof. We have Taylor expansion (using the Lagrange form for the remainder) hp(¯ µ(M) i) = hp(µi) + h0p(µi)¯ µ(M) i−µi+1 2! h00 p(m)¯ µ(M) i−µi2, for some mbetween ¯ µ(M) iand µi. This provides M1/2 hp(¯ µ(M) i)−hp(µi) h0p(µi)Var(b µi)1/2 =M1/2 ¯ µ(M) i−µi Var(b µi)1/2 +1 2!h0p(µi)Var(b µi)1/2 M1/2 h00 p(m)¯ µ(M) i−µi2. The first term on the right-hand side converges in distribution to the standard Gaussian distribution as M→∞ , because of the CLT Equation (11) . Therefore, the claim follows by proving that the last term converges in probability to zero. Consider the event AM={|¯ µ(M) i−µi|>M−3/8}= Risks 2020,8, 83 15 of 26 from the EDF may be more appropriate. In any of these other model choices the proposed methodology works similarly, only minor changes to the Rcode in Listing A1 will be necessary: • Firstly, the complexity of the network architecture should be adapted to the problem, both the sizes of the hidden layers and the depth of the network will vary with the complexity of the problem and the complexity of the regression function. In general, the bigger the network the better the approximation capacity (this follows from the universality theorems). This says that the chosen network should have a certain complexity otherwise it will be not sufficiently flexible to approximate the true (but unknown) regression function. On the other hand, for computational reasons, the chosen network should not be too large. • Secondly, different distributions will require different choices of loss functions on line 25 of Listing A1. Some loss functions are already implemented in the keras library, others will require custom loss functions. An example of a custom loss function implementation is given in Listing 2 of Delong et al. (2020). In general, the loss function of any density that allows for an EDF representation can be implemented in keras. • Thirdly, we could add more specialized layers to the network architecture of Listing A1. For instance, we could add dropout layers after each hidden layer on lines 15–17 of Listing A1, for dropout layers see Srivastava et al. (2014). Dropout layers add an additional element of randomness during training, because certain network connections are switched off (at random) for certain training steps when using dropout. This switching off of connections acts as regularization during training because it prevents certain neurons from over-fitting to special tasks. In fact, under certain assumptions one can prove that dropout acts similarly to ridge regularization, see Section 18.6 in Efron and Hastie (2016). In our study we refrain from using dropout. 5.5. Comparison of Different Network Calibrations The issue with the network result on line (d) of Table 2now is that it involves quite some randomness: (R1) we randomly split learning data Dinto training data D(−)and validation data V; (R2) we randomly split training data D(−) into mini-batches of size 5000 (to more efficiently calculate gradient descent steps); and (R3) we randomly choose the starting point of the gradient descent algorithm; the default initialization in keras is the glorot_uniform initialization which involves simulation from uniform distributions. Changing seeds in these points (R1)–(R3) will change the network regression parameter estimate b β(j) and, hence, the results. For this study we explore the randomness of changing the seeds in (R1)–(R3). However, one may also be interested in sensitivities of results when changing sizes of the mini-batches, hyper-parameters of the gradient descent methods, or when using stratified versions of validation data. We will not explore this here. We run the above calibration procedure under identical choices of all hyper-parameters, but we choose different seeds for the random choices (R1)–(R3). The boxplot in Figure 2shows the in-samples losses L(D ; b β(j)) and out-of-samples losses L(T ; b β(j)) over 400 network calibrations b β(j) by only randomly changing the seeds in the above mentioned points (R1)–(R3). We note that these losses have a rather large range which indicates that results of single network calibrations are not very robust. We can calculate empirical mean and standard deviation for the 400 seeds jgiven by L(T;b β(1:400)) = 1 400 400 ∑ j=1L(T;b β(j))and v u u t1 399 400 ∑ j=1L(T;b β(1:400))−L(T;b β(j))2 . The first gives an estimate for the expected generalization loss Equation (8) averaged over the corresponding portfolios. We emphasize in notation b β(1:400) that we do not average over network Risks 2020,8, 83 16 of 26 parameters, but over deviance losses on individual network parameters b β(j) . The resulting numbers are given on line (e) of Table 2. This shows that the early stopped network calibrations have quite significant differences, which motivates the study of the nagging predictor. ● ● ● ● ● ● ● ● ● ● 30.1 30.2 30.3 30.4 30.5 in−sample: boxplot over 400 calibrations in−sample losses 400 calibrations selected calibration ● ● ● ● ● ● ● 31.4 31.5 31.6 31.7 out−of−sample: boxplot over 400 calibrations out−of−sample losses 400 calibrations selected calibration Figure 2. In-sample losses ( top ) and out-of-sample losses ( bottom ) over 400 different runs of the gradient descent algorithm with different seeds j= 1, . . . , 400, the blue and red lines show the network regression model on line (d) of Table 2with seed j=1. Figure 3gives a scatter plot of in-sample and out-of-sample losses over the 400 different runs of the gradient descent fitting (complemented with a natural cubic spline). We note in this example that both small and big in-sample losses do not lead to favorable generalization, small in-sample losses may indicate over-fitting and big out-of-sample losses may indicate the that early stopping rule has not found a good model based on the validation data V. Risks 2020,8, 83 17 of 26 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● 30.1 30.2 30.3 30.4 30.5 31.4 31.5 31.6 31.7 scatter plot of in−sample and out−of−sample losses in−sample losses out−of−sample losses ●400 calibrations cubic spline Figure 3. Scatter plot of in-sample losses versus out-of-sample losses over 400 different runs of the gradient descent algorithm with different seeds j= 1, . . . , 400, the orange line provides a cubic spline approximation. 5.6. Nagging Predictor We are now in the situation where we are able to calculate the nagging predictors ¯ ¯ µ(M) t over the test data set T . For M→∞ this provides us with empirical counterparts of Propositions 3and 4. We therefore consider for M≥1 the sequence of out-of-sample losses, see Equation (12), L(T;¯ ¯ µ(M) t=1,...,m) = 1 m m ∑ t=1 δ(Y† t,¯ ¯ µ(M) t), where ¯ ¯ µ(M) t are the nagging predictors Equation (13) received from the i.i.d. sequence b µ(j) t=µ(x† t , b β(j)) and, again, where the i.i.d. property is conditional on D and refers to the i.i.d. starting points chosen to start the gradient descent fitting algorithm. Figure 4gives the out-of-sample losses of the nagging predictors L(T ; ¯ ¯ µ(M) t=1,...,m) for M= 1, . . . , 100. Most noticeable is that nagging leads to a substantial improvement in out-of-sample losses, for M→∞ the out-of-sample loss converges to 31.272 which is much smaller then the figures reported on lines (b)–(e) of Table 2, in fact, the out-of-sample loss decreases by more than 3 standard deviations (as reported on line (e) of Table 2). From this we conclude that nagging helps to improve the predictive model substantially. From Figure 4we also observe that this convergence mainly takes place over the first 20 aggregating steps in our example. That is, we need to aggregate over roughly 20 network calibrations b β(j) to get the maximal predictive power. The dotted orange lines in Figure 4give corresponding 1 standard deviation confidence bounds (received by repeating the nagging procedure for different seeds). To get sufficiently small confidence bounds in our example we need to average over roughly 40 network calibrations. Risks 2020,8, 83 18 of 26 0 20 40 60 80 100 31.30 31.35 31.40 31.45 31.50 31.55 nagging predictors for M>=1 index M out−of−sample losses nagging predictor 1 standard deviation Figure 4. Out-of-sample losses of nagging predictors ¯ ¯ µ(M) taveraged over Tfor M=1, . . . , 100. Conclusion 1. Nagging helps to sufficiently improve out-of-sample performance of network regression models. In our example we need to average over 20 different calibrations to get optimal predictive power, and we need an average of 40 calibrations to ensure that we cannot further improve this predictive power. 5.7. Pricing of Individual Insurance Policies In the previous sections we have proved that nagging successfully improves models. However, all these considerations have mainly been based on portfolio considerations, i.e., our focus has been on the question whether the insurance company has a model that gives good predictive power to predict the portfolio claim amount. For insurance pricing, this is not sufficient because we should also ensure that we have robustness of prices on an individual insurance policy level. In this section we analyze by how much individual insurance policy prices may differ if we select two different network calibrations b β(j) and b β(j0) . This will tell us whether aggregating over 20 or 40 network calibrations is sufficient as stated in Conclusion 1. Naturally, we expect that we need to average over more networks because the former statement includes an average over T , i.e., over m= 67,801 insurance policies (though there is dependence between these policies because they simultaneously use the same network parameter estimate b β(j)). To analyze this question on individual insurance policies we calculate for each policy t= 1, . . . , m of the test data T the nagging predictor ¯ ¯ µ(M) t over M= 400 different network calibrations b β(j) , j= 1, . . . , M , and we calculate the empirical coefficients of variation in the individual network predictors given by d CoVt=bσt ¯ ¯ µ(M) t =r1 M−1∑M j=1b µ(j) t−¯ ¯ µ(M) t2 ¯ ¯ µ(M) t , (15) these are the empirical standard deviations normalized by the corresponding mean estimates. We plot these coefficients of variations d CoVt in Figure 5(lhs) against the nagging predictors ¯ ¯ µ(M) t for each single insurance policy t= 1, . . . , m (out-of-sample on T ). Figure 5(rhs) shows the resulting histogram. We observe that on most insurance policies (73%) we have a coefficient of variation of less than 0.2, however, on 11 of the m= 67,801 insurance policies we have a coefficient of variation bigger then 1. Thus, for the latter, if we average over 400 different network calibrations b β(j) we still have an uncertainty of 1 /√400 = 0.05, i.e., the prices have a precision of 5% to 10% in these latter cases (this is always conditional given D ). From this we conclude that on individual insurance policies we need to aggregate over a considerable number of networks to receive stable network regression prices. Risks 2020,8, 83 19 of 26 histogram of coefficients of variations coefficient of variations frequency 0.1 0.2 0.3 0.4 0 2000 4000 6000 8000 10000 Figure 5. Coefficients of variation d CoVt on individual policy predictions t= 1, . . . , m : (lhs) scatter plot against nagging predictors ¯ ¯ µ(M) t, and (rhs) histogram. We show these 11 insurance policies with a coefficient of variation bigger 1 in Listing 2. Striking is that all these insurance policies have vehicle age VehAge = 0. In Figure 6(top row) we show a scatter plot and histogram of these insurance policies, and Figure 6(bottom row) gives the corresponding graphs for VehAge > 0. We indeed confirm that mainly the policies with VehAge = 0 are difficult to price. From the histogram in Figure 5(rhs) we receive in total 1423 policies that have a coefficient of variation d CoVt bigger than 0.4, and in view of Figure 6(rhs), 1163 of these insurance policies have VehAge =0. Listing 2. Policies with high coefficients of variation d CoVt. 1Area VehPower VehAge DrivAge BonusMalus VehBrand VehGas Density Region 2A 6 0 51 50 B3 Diesel 2.71 R21 3A 6 0 51 50 B3 Diesel 2.71 R21 4B 9 0 30 125 B3 Regular 4.32 R26 5E 15 0 75 67 B14 Regular 8.38 R72 6B 6 0 29 60 B3 Diesel 4.30 R21 7A 10 0 29 60 B13 Regular 2.08 R24 8E 9 0 31 125 B4 Diesel 8.35 R11 9A 7 0 69 50 B14 Diesel 3.83 R82 10 A 10 0 59 50 B1 Diesel 3.33 R21 11 A 10 0 59 50 B1 Diesel 3.33 R21 12 A 10 0 59 50 B1 Diesel 3.33 R21 To better understand this uncertainty in VehAge = 0, we provide some empirical plots for those vehicles with age 0; these plots are taken from Noll et al. (2018). Figure 7shows the total exposure per vehicle age (top-lhs), the empirical frequency per vehicle age (top-rhs), the vehicle ages per vehicle brands (bottom-lhs) and total exposures per vehicle brand (bottom-rhs). From these graphs we observe that vehicle age 0 has a completely different frequency than all other vehicle ages (Figure 7, top-rhs). Moreover, vehicle age 0 is dominated by vehicle brand B12 (Figure 7, bottom). It seems in view of Listing 2that this configuration results in quite some uncertainty in frequency prediction, note that 1011 of the 1163 insurance policies with vehicle age 0 and a coefficient of variation bigger than 0.4 are cars of vehicle brand B12 . Naturally, in a next step we would have to analyze vehicle age 0 and vehicle brand B12 in more depth, we suspect that these could likely be rental cars (or some other special cases). Unfortunately, there is no further information available for this data set that allows us to investigate this problem more thoroughly. We remark that König and Loser (2020) come to a similar conclusion for cars of vehicle brand B12 with vehicle age 0, see regression tree in Chapter 1.2 of König and Loser (2020) where furthermore the split w.r.t. vehicle gas is important. Risks 2020,8, 83 20 of 26 histogram of coefficients of variations (only VehAge=0) coefficient of variations frequency 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0 200 400 600 800 1000 1200 histogram of coefficients of variations (VehAge>0) coefficient of variations frequency 0.1 0.2 0.3 0.4 0 2000 4000 6000 8000 10000 Figure 6. ( Top row ) Only vehicle age VehAge = 0, ( bottom row ) vehicle age VehAge > 0: coefficients of variation d CoVt on individual policy predictions t= 1, . . . , m : (lhs) scatter plot against nagging predictors ¯ ¯ µ(M) t, and (rhs) histogram. 5.8. Meta Network Regression Model From the previous sections we conclude that the nagging predictor substantially improves the predictive model. Nevertheless, the nagging predictor may not be fully satisfactory in practice. The difficulty is that it involves aggregating over M= 400 predictors b µ(j) i for each policy i . Such blended predictors (and models) are not easy to maintain nor is it simple to study model properties, updating models with new observations, etc., for instance, if we have a new insurance policy having covariate x , then we need to run this new insurance policy through all M= 400 networks to receive the price. For this reason we propose to build a meta model that fits a new network to the nagging predictors ¯ ¯ µ(M) i , i= 1, . . . , n . In the wider machine learning literature, building such meta models is also referred to as “model distillation”, see Hinton et al. (2015). Since these nagging predictors are aggregated predictors over M network models, and since the network regression functions themselves are smooth functions in input variables (at least in the continuous features), the nagging predictors describe smooth surfaces. Therefore, it is comparably simple to fit a network to the smooth surface described by nagging predictors ¯ ¯ µ(M) i , i= 1, . . . , n , and over-fitting will not be an issue. To build this meta model we use exactly the same network architecture as above, and as outlined in Listing A1 Risks 2020,8, 83 21 of 26 in the Appendix A. The only parts that we change are the loss function and the response variables. We replace the original claim count responses Yi by the nagging predictors ¯ ¯ µ(M) i , and for the loss function we choose the square loss function. We can either choose an unweighted square loss function or we can weight the individual observations with the inverse standard deviation estimates 1 /bσi , see Equation (15) , the latter reflects the idea of giving more weight to predictors that have less uncertainty. In fact, the latter is a little bit similar to Zhou et al. (2002) where the authors argue that equally weighting over network predictors is non-optimal. Our approach is slightly different here because we do not weight networks, but we weight nagging predictors of individual insurance policies to receive better rates of convergence for meta model training. 0 1 2 3 4 5 6 7 8 9 11 13 15 17 19 total volumes per vehicle age groups vehicle age groups exposure 0 5000 10000 15000 20000 25000 30000 ● ●●●●●●●●●●●●●●●●●●●● 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 observed frequency per vehicle age groups vehicle age groups frequency ● ●●●●●●● ●●●●● ●●●●● ●●● 0 1 2 3 4 5 6 7 8 9 11 13 15 17 19 B1 B10 B11 B12 B13 B14 B2 B3 B4 B5 B6 10 9 8 7 6 5 4 3 2 1 0 vehicle age among vehicle brands car brand groups relative frequency 0.0 0.2 0.4 0.6 0.8 1.0 B1 B10 B11 B12 B13 B14 B2 B3 B4 B5 B6 total volumes per car brand groups car brand groups exposure 0 20000 40000 60000 80000 Figure 7. Empirical plots of VehAge s: ( top lhs ) exposure per vehicle ages, ( top rhs ) marginal frequencies per vehicle ages, ( bottom lhs ) interaction vehicle age and vehicle brand, ( bottom rhs ) exposures per vehicle brand. We present the gradient descent fitting performance in Figure 8and the resulting in-sample and out-of-sample losses in Table 3. We conclude that the weighted version of the square loss has better convergence properties in gradient descent fitting, and the resulting model has a better loss performance, compare lines (g1) and (g2) of Table 3. The resulting meta model on line (g2) is slightly worse than the nagging predictor model on line (f), however, substantially better than the individual Risks 2020,8, 83 22 of 26 network models (d)–(e) and much more easy in handling than the nagging predictor. For this reason, we are quite satisfied by the meta model, and we propose to hold on to this model for further analysis and insurance pricing. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 200 400 600 800 1000 0.0000 0.0005 0.0010 0.0015 0.0020 stochastic gradient descent algorithm training epochs deviance loss ● training loss validation loss ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0 200 400 600 800 1000 0.00 0.01 0.02 0.03 0.04 0.05 0.06 stochastic gradient descent algorithm training epochs deviance loss ● training loss validation loss Figure 8. Gradient descent fitting of meta model: (lhs) unweighted square loss function, (rhs) weighted square loss function. Table 3. (d)–(f) Network regression models from Table 2compared to meta models (g1)–(g2); losses are in 10−2. In-Sample Out-of-Sample Loss on DLoss on T (d) network regression model (seed j=1) 30.184 31.464 (e) average over 400 network calibrations 30.230 (0.089) 31.480 (0.061) (f) nagging predictor for M=400 30.060 31.272 (g1) meta network model (un-weighted) 30.260 31.342 (g2) meta network model (weighted) 30.257 31.332 Remark 4. We have mentioned in Remark 3that model implementation will slightly change if we move from our Poisson example to another distributional example, for instance, the loss function has to be adapted to the choice of the distribution function. The last step of fitting a meta model, however, does not require any changes. Note that the meta model fits a regression function to another regression function (which is the nagging predictor in our case). This fitting uses the square loss as the objective function, and the nature of the distribution function of the observations is not important because we fit a deterministic function to another deterministic function here. In Figure 9we present the two predictors in a scatter plot, we observe that they are reasonably equal, the biggest differences are highlighted in blue color, and they refer to policies with vehicle age 0, that is, the feature component within the data that is the most difficult to fit with the network model. Risks 2020,8, 83 23 of 26 Figure 9. Comparison of nagging predictors and meta model predictors. We provide a concluding analysis of this meta model. As emphasized in Section 7 of Richman et al. (2019) , we analyze and try to understand the representations learned in the last hidden layer of the network over all insurance policies. In view of Equation (1) this means that we study the learned representations for t=1, . . . , mgiven by x† t7→ ˘x† t=z(d)◦···◦z(1)(x† t)∈Rqd+1, to which the GLM step with readout parameter β(d+1)∈Rqd+1 is applied to. In our example of Listing A1 we have d= 3 hidden layers with q3= 10 hidden neurons in this last hidden layer. In order to illustrate the learned representations we apply a principal component analysis (PCA) to these learned representations (˘xt)t=1,...,m . This provides us with the principal component (PC) scores for each policy t and the corresponding singular values λq3≥. . . ≥λ1≥ 0. Note that these PC scores are ordered w.r.t. importance measured by the singular values. We scale the individual PC scores k= 1, . . . , q3 with √λk to reflect this importance. In Figure 10 (top row) we plot these scaled PC scores for the first 4 principal components and averaged per vehicle age (lhs), driver age (middle), and bonus malus level (rhs). We remark that these 3 feature components are the most predictive ones according to the descriptive analysis in Noll et al. (2018). From Figure 10 (top row) we observe that the first principal component (in red color) reflects the marginal empirical frequencies illustrated in Figure 10 (bottom row), thus, the first PC reflects the influence of each of these three explanatory variables on the regression function, and higher order PCs are used to refine this picture and model interactions between feature components. Risks 2020,8, 83 24 of 26 0 5 10 15 20 −0.10 −0.05 0.00 0.05 0.10 0.15 0.20 0.25 PC scores for vehicle age vehicle age scaled PC score PC score 1 PC score 2 PC score 3 PC score 4 20 30 40 50 60 70 80 90 −0.10 −0.05 0.00 0.05 0.10 0.15 0.20 0.25 PC scores for driver age driver age scaled PC score PC score 1 PC score 2 PC score 3 PC score 4 60 80 100 120 140 −0.10 −0.05 0.00 0.05 0.10 0.15 0.20 0.25 PC scores for bonus malus level bonus malus level scaled PC score PC score 1 PC score 2 PC score 3 PC score 4 ● ●●●●●●●●●●●●●●●●●●●● 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 observed frequency per vehicle age groups vehicle age groups frequency ● ●●●●●●● ●●●●● ●●●●● ●●● 0 1 2 3 4 5 6 7 8 9 11 13 15 17 19 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 observed frequency per driver's age groups driver's age groups frequency 18 23 28 33 38 43 48 53 58 63 68 73 78 83 88 ● ●●●● ● ● ● ● ● ● 0.0 0.1 0.2 0.3 0.4 0.5 0.6 observed frequency per bonus−malus level groups bonus−malus level groups frequency 50 60 70 80 90 100 110 120 130 140 150 Figure 10. ( top ) Scaled principal component (PC) scores on the first 4 principal components averaged over the corresponding vehicle ages (lhs), driver ages (middle), and bonus males levels (rhs); ( bottom ) observed (empirical) marginal frequencies of vehicle ages (lhs), driver ages (middle), and bonus males levels (rhs). 6. Conclusions and Outlook This work has defined the nagging predictor which produces accurate and stable portfolio predictions on the basis of random network calibrations, and it provided convergence results in the context of Tweedie’s compound Poisson generalized linear models. Focusing on an example in motor third-party liability insurance pricing, we have shown that stable portfolio results are achieved after 20 network training runs, and by increasing the number of network training runs to 400, we have shown that quite stable results are produced at the level of individual policies, which is an important requirement for the use of networks for insurance pricing, and more general actuarial tasks. The coefficient of variation of the nagging predictor is shown to be a useful data-driven metric for measuring the relative difficulty with which a network is able to fit to individual training examples, and we have used it to calibrate an accurate meta network which approximates the nagging predictor. Whereas this work has examined the stability of network predictions in the case of a portfolio at a point in time, another important aspect of consistency within insurance is stable pricing over time, thus, future work could consider methods for stabilizing network predictions as new information becomes available. Author Contributions: Both authors (R.R., M.V.W.) have equally contributed to the this project, this concerns the concept, the methodology, the writing, and the numerical analysis of this project. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Conflicts of Interest: The authors declare no conflict of interest. Appendix A. RCode Listing A1 provides the Rcode used to fit the data.