Log-normal or over-dispersed poisson?

Harnau, Jonas Article Log-normal or over-dispersed poisson? Risks Provided in Cooperation with: MDPI – Multidisciplinary Digital Publishing Institute, Basel Suggested Citation: Harnau, Jonas (2018) : Log-normal or over-dispersed poisson?, Risks, ISSN 2227-9091, MDPI, Basel, Vol. 6, Iss. 3, pp. 1-37, https://doi.org/10.3390/risks6030070 This Version is available at: https://hdl.handle.net/10419/195862 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 Log-Normal or Over-Dispersed Poisson? Jonas Harnau ID Department of Economics, University of Oxford & Oriel College, Oxford OX1 4EW, UK; [email protected] Received:18 June 2018; Accepted: 6 July 2018; Published: 9 July 2018 Abstract: Although both over-dispersed Poisson and log-normal chain-ladder models are popular in claim reserving, it is not obvious when to choose which model. Yet, the two models are obviously different. While the over-dispersed Poisson model imposes the variance to mean ratio to be common across the array, the log-normal model assumes the same for the standard deviation to mean ratio. Leveraging this insight, we propose a test that has the power to distinguish between the two models. The theory is asymptotic, but it does not build on a large size of the array and, instead, makes use of information accumulating within the cells. The test has a non-standard asymptotic distribution; however, saddle point approximations are available. We show in a simulation study that these approximations are accurate and that the test performs well in finite samples and has high power. Keywords: non-nested testing; encompassing; chain-ladder 1. Introduction Which is the better chain-ladder model for claim reserving: over-dispersed Poisson or log-normal? While the expert may have a go-to model, the answer should be informed by the data. Choosing the wrong model could substantially influence the quality of the reserve forecast. Yet, so far, no statistical theory is available that supports the actuary in his/her decision and that allows him/her to make a solid argument in favour of either model. We develop a test that can distinguish between over-dispersed Poisson and log-normal data generating processes, both of which have a long history in claim reserving. The test exploits that the former model fixes the variance to mean ratio across the array, while the latter assumes a common standard deviation to mean ratio. Consequently, the test statistic is based on estimators for the variation in the respective models. The idea is drawn from the econometric literature on encompassing. Intuitively, the test asks whether the null-model can accurately predict the behaviour of the rival model’s variation estimator when the null-model is true. The over-dispersed Poisson model is appealing since it naturally pairs with Poisson quasi-likelihood estimation, replicating the popular chain-ladder technique in run-off triangles (Kremer 1985, pp. 130) . Furthermore, this model makes for an appealing story due to its relation to compound Poisson distributions. Such distributions give the aggregate incremental claims an interpretation as the sum over a Poisson number of claims with random individual claim amounts (Beard et al. 1984, Section 3.2). A popular method to generate distribution forecasts for the over-dispersed Poisson model is bootstrapping (England and Verrall 1999;England 2002). While in widespread use, there is so far no theory proving the validity of the bootstrap in this setting. Furthermore, in some settings, the method seems to produce unsatisfactory results. Recently, Harnau and Nielsen (2017) developed a theory that gives the over-dispersed Poisson model a rigorous statistical footing. They propose an asymptotic framework based on infinitely-divisible distributions that keeps the dimension of the data array fixed and instead builds on large cell means. This resolves the incidental parameter problem (Lancaster 2000;Neyman and Scott 1948) that renders a standard asymptotic theory based on a large array invalid and arises since the number of parameters grows with the size of the array. The class of infinitely-divisible distributions includes compound Poisson Risks 2018,6, 70; doi:10.3390/risks6030070 www.mdpi.com/journal/risks Risks 2018,6, 70 2 of 37 distributions, which are appealing in an insurance context, as noted above. We can then interpret large cell means as the result of a large latent underlying number of claims. Other infinitely-divisible distributions that can be reconciled with the over-dispersed Poisson structure include Poisson, gamma and negative binomial. The intuition for the theory by Harnau and Nielsen (2017) is that the array is roughly normally distributed for large cell means, so that the results remind us of a classical analysis of variance (ANOVA) setting. Harnau and Nielsen (2017) show that Poisson quasi-likelihood estimators are t -distributed, and F -tests based on Poisson likelihoods can be used to test for model reduction, such as for the absence of a calendar effect. Finally, chain-ladder forecast errors are t -distributed, giving rise to closed-form distribution forecasts including for aggregates, such as the reserve or cash-flow. In their simulations, Harnau and Nielsen (2017) find that while the bootstrap (England and Verrall 1999;England 2002) matches the true forecast error distribution better on average, the t -forecast produces fewer outliers and appears more robust. Building on the asymptotic framework put forward by Harnau and Nielsen (2017), Harnau (2018a) proposes misspecification tests for two crucial assumptions of the over-dispersed Poisson model. First, the variance to mean ratio is assumed to be common across the array. Second, accident effects are not allowed to vary over development years and vice versa. To check for a violation of these assumptions, Harnau (2018a) suggests splitting the run-off triangle into sub-samples and then testing whether a reduction from individual models for each sub-sample to a single model for the full array can be justified. While the idea of splitting the sample is borrowed from time-series econometrics (Chow 1960) , the theory for a reduction to a single model is again reminiscent of an ANOVA setting. A classical Bartlett test (Bartlett 1937) can be used to assess whether we can justify common variance to mean ratios. This is followed by an independent F -test for the absence of breaks in accident and development effects. Again, the asymptotics needed to arrive at these results keep the dimension of the array fixed, growing instead the cell means. Harnau (2018a) also shows that these misspecification tests can be used in a similar fashion in a finite sample log-normal model. The log-normal model introduced by Kremer (1982), who relates it to the ANOVA literature, features a predictor structure that is reminiscent of the classical chain-ladder. Verrall (1994) refers to this as the chain-ladder linear model, while Kuang et al. (2015) use the term geometric chain-ladder. The latter authors show that the maximum likelihood estimators in the log-normal model can be interpreted as development factors of geometric averages, compared to an interpretation of arithmetic averages arising for the classical chain-ladder. An advantage of the log-normal model is that an exact Gaussian distribution theory applies to the maximum likelihood estimators. However, since these estimators are computed on the log scale, a bias is introduced on the original scale. Verrall (1991) tackles this issue and derives unbiased estimators for the mean and standard deviation on the original scale. One issue for full distribution forecasts in the log-normal model is that the insurer is usually not interested in forecasts for individual cells, but rather for cell sums such as the reserve or the cash-flow. However, the log-normal distribution is not closed under convolution, so that cell sums are not log-normally distributed. Recently, Kuang and Nielsen (2018) proposed a theory that includes closed-form distribution forecasts for cell sums, such as the reserve, in the log-normal model, thus remedying one of its drawbacks. Kuang and Nielsen (2018) combined the insight by Thorin (1977) that the log-normal distribution is infinitely divisible and the asymptotic framework by Harnau and Nielsen (2017). Based on this, they propose a theory for generalized log-normal models, a class that nests the log-normal model, but is not limited to it. In particular, the distribution is not assumed to be exactly log-normal, but merely needs to be infinitely divisible with a moment structure close to that of the log-normal model. The asymptotics in this framework again leave the dimension of the array untouched to avoid an incidental parameter problem. In contrast to the theory for large cell means in the over-dispersed Poisson model, results are now for small standard deviation to mean ratios. Risks 2018,6, 70 3 of 37 For the generalized log-normal model, Kuang and Nielsen (2018) show that least squares estimators computed on the log scale are asymptotically t -distributed, and simple F -tests based on the residual sum of squares can be used to test for model reduction. Reassuringly, these results match the exact results in a log-normal model. Beyond that, they also prove that forecast errors on the original scale are asymptotically t -distributed so that distribution forecasting for cell sums is straightforward. Further, they show that the misspecification tests by Harnau (2018a) are asymptotically valid for the generalized log-normal model, just as they were in finite samples for the log-normal model. We remark that besides over-dispersed Poisson and log-normal models, there exist a number of reserving models that we do not consider further in this paper. England and Verrall (2002) give an excellent overview. Perhaps the most popular contender is the “distribution-free” model by Mack (1993) . This model also replicates the classical chain-ladder, but differs from the over-dispersed Poisson model. Mack (1993) derives the expression for forecast standard errors. However, so far, no full distribution theory exists for this model. With a range of theoretical results in place for over-dispersed Poisson and (generalized) log-normal models, discussed further in Section 3, a natural question is when we should employ which model. The misspecification tests by Harnau (2018a) seem like a natural starting point. For example, if we can reject the specification of the log-normal, but not the over-dispersed Poisson model, the latter seems preferable. However, the misspecification tests may not always have enough power to make this distinction, as we show in Section 2. Since generalized log-normal and over-dispersed Poisson models are not nested, a direct test between them is not trivial. Cox (1961,1962) introduced a theory for non-nested hypothesis testing with a null model. Vuong (1989) provided a theory for non-nested model selection without a null model; in selection, the goal is to choose the better, not necessarily the true, model. However, both procedures are likelihood based, so that the results are not applicable here, since we did not specify exact distributions and, thus, do not have likelihoods available. Given the lack of likelihoods for the models, we look to the econometric encompassing literature for inspiration. The theory for encompassing allows for a more general way of non-nested testing. As Mizon and Richard (1986) put it, “Among other criteria, it seems natural to ask whether a specific model, say M1 , can mimic the DGP [data generating process], in that statistics which are relevant within the context of another model, M2 say, behave as they should were M1 the DGP.” The encompassing literature originates from Hendry and Richard (1982) and Mizon and Richard (1986); for a less technical introduction, see Hendry and Nielsen (2007, Section 11.5). Ermini and Hendry (2008) applied the encompassing principle in a time-series application. They tested whether disposable income is better modelled on the original scale or in logs. Taking the log model as the null hypothesis, they evaluated whether the log model can predict the behaviour of estimators for the mean and variance of the model on the original scale. Building on the encompassing literature, we find the distribution of the over-dispersed Poisson model estimators under a generalized log-normal data generating process and vice versa. It turns out that both Poisson quasi-likelihood and log data least squares estimators for accident and development effects are asymptotically normal, regardless of the data generating process. Differences arise in the second moments. This manifests in the limiting distributions of the variation estimators. While these are asymptotically χ2 under the correct model, their distribution is a non-standard quadratic form of normals under the rival model. However, these distributions involve the unknown dispersion parameter, which needs to be estimated. Employing the variation estimator of the correct model for this purpose, we arrive at a test statistic with a non-standard asymptotic distribution: the ratio of dependent quadratic forms. Saddle point approximations to such distributions are available (Butler and Paolella 2008;Lieberman 1994). Further, we can show that the power of the tests originates from variation in the means across cells. This is intuitive given that the main difference between the models disappears when all means are identical; then, both standard deviation to the mean and variance to mean ratios are constant across the array. These findings are collected in Section 4. Risks 2018,6, 70 4 of 37 With the theoretical results for encompassing tests between over-dispersed Poisson and generalized log-normal models in place, we show that they perform well in a simulation study. First, we demonstrate that saddle point approximations to the limiting distributions of the statistics work very well. Second, we tackle an issue that disappears in the limit: we have the choice between a number of asymptotically-identical estimators that generally differ in finite samples. Simulations reveal substantial heterogeneity in finite sample performance, but also show that some choices generally do well. Third, we show that the tests have high power for parameterizations we may realistically encounter in practice. We also find that power grows quickly with the variation in the means. The simulation study is in Section 5. Having convinced ourselves that the tests do well in simulations, we demonstrate their application in a range of empirical applications in Section 6. First, we revisit the empirical illustration of the problem from the beginning of the paper. We show that the test has no problem rejecting one of the two rival models. Second, we consider an example that perhaps somewhat cautions against starting with a model that may be misspecified to begin with. In this application, dropping a clearly needed calendar effect turns the results of the encompassing tests upside down. Third, taking these insights into account, we implement a testing procedure that makes use of a whole range of recent results: deciding between the over-dispersed Poisson and generalized log-normal model, evaluating misspecification and testing for the need for a calendar effect. We conclude the paper with a discussion of potential avenues for future research in Section 7. These include further misspecification tests, a theory for the bootstrap and empirical studies assessing the usefulness of the recent theoretical developments in applications. 2. Empirical Illustration of the Problem We illustrate in an empirical example that the choice between the over-dispersed Poisson and (generalized) log-normal model is not always obvious. Table 1shows a run-off triangle taken from Verrall et al. (2010, Table 1) with accident years i in the rows and development years j in the columns. Calendar years k=i+j−1 are on the diagonals. Table 1. Run-off triangle taken from Verrall et al. (2010) with an indication for splitting into sub-samples corresponding to the first and last five accident years. Accident years i in the rows; development years jin the columns. i, j 1 2 3 4 5 6 7 8 9 10 1 451,288 339,519 333,371 144,988 93,243 45,511 25,217 20,406 31,482 1729 2 448,627 512,882 168,467 130,674 56,044 33,397 56,071 26,522 14,346 - 3 693,574 497,737 202,272 120,753 125,046 37,154 27,608 17,864 - - 4 652,043 546,406 244,474 200,896 106,802 106,753 63,688 - - - 5 566,082 503,970 217,838 145,181 165,519 91,313 - - - - 6 606,606 562,543 227,374 153,551 132,743 - - - - - 7 536,976 472,525 154,205 150,564 - - - - - - 8 554,833 590,880 300,964 - - - - - - - 9 537,238 701,111 - - - - - - - - 10 684,944 - - - - - - - - - While Kuang et al. (2015) and Harnau (2018a) model the data in Table 1as log-normal, it is not obvious whether a log-normal or an over-dispersed Poisson model is more appropriate. In a log-normal model, the aggregate incremental claims Yij are independent: MLN : log(Yij) = N(αi+βj+δ,ω2) Risks 2018,6, 70 5 of 37 where α and β are accident and development effects, respectively. On the original scale, this implies that: MLN =⇒E(Yij) = exp αi+βj+δ+ω2 2,sd(Yij) E(Yij)=qexp(ω2)−1. (1) Thus, the standard deviation to mean ratio, as well as the log data variance are common across cells. If we instead chose an over-dispersed Poisson model, we would maintain the independence assumption and specify the first two moments of the claims Yij as: MODP :E(Yij) = exp(αi+βj+δ),var(Yij) E(Yij)=σ2. Thus, the variance to mean ratio σ2is identical for all cells. To choose between the two models, we could take the misspecification tests by Harnau (2018a) as a starting point. To implement the tests, the data are first split into sub-samples. For the Verrall et al. (2010) data, Harnau (2018a) considers a split into two sub-samples consisting of cells relating to the first and last five accident years, as illustrated in Table 1. The idea is then to test for common parameters across the sub-samples. In the log-normal model, we first perform a Bartlett test for common log data variances ω2 across sub-samples and, if this is not rejected, an F -test for common accident and development effects. Similarly, in the over-dispersed Poisson model, we first test for common over-dispersion σ2and then again for common accident and development effects. If one of the models is flagged as misspecified, but not the other, the choice becomes obvious. However, in this application, we cannot reject either model based on these tests. For the log-normal model, the Bartlett test for common log data variances yields a p-value of 0.09 and the F -test for common effects a p-value of 0.91; the p-values for the equivalent tests in the over-dispersed Poisson model are 0.78 and 0.64, respectively. Therefore, the question remains: Which model should we choose? 3. Overview of the Rival Models We first discuss two common elements of the rival models, namely the data structure and the chain-ladder predictor and its identification. Then, we in turn state assumptions, estimation and known theoretical results for the over-dispersed Poisson and the generalized log-normal chain-ladder model. 3.1. Data We assume that we have data for a run-off triangle of aggregate incremental claims. We denote the claims for accident year i and development year j by Yij . Further, we count calendar years with an offset, so the calendar year k=i+j− 1. Then, we can define the index set for a run-off triangle with I accident, development and calendar years by: I={(i,j): 1 ≤i,j,k≤I}. We define the number of observations in I as n . We could also allow for data in a generalized trapezoid as defined by Kuang et al. (2008) without changing the results of the paper. Loosely, generalized trapezoids allow for an unbalanced number of accident and development years, as well as missing calendar years both in the past and the future. 3.2. Identification We briefly discuss the identification problem of the chain-ladder predictor αi+βj+δ that is common to both over-dispersed and generalized log-normal models. Kremer (1985) showed that based on this predictor, Poisson quasi-likelihood estimation replicates the classical chain-ladder point forecasts in a run-off triangle. Risks 2018,6, 70 6 of 37 The identification problem is that for any aand b, µij =αi+βj+δ= (αi+a) + (βj+b) + (δ−a−b) where αi and βj are accident and development effects, respectively. Thus, no individual effect is identified. Several ad hoc identification methods are available; for example, we could set ∑iαi=∑jβj=0 .Kuang et al. (2008) suggest a parametrization that is canonical in a Poisson model and allows for easy counting of degrees of freedom. The idea is to re-write the linear predictor in terms of a level and deviations from said level as: µij =µ11 + I ∑ s=2 1(i≤s)∆αs+ I ∑ t=s 1(j≤s)∆βs. Thus, we can write: µij =x0 ijξ where the design xij and identified parameter ξare given by: xij = (1, 1(i≤2), . . . , 1(i≤I), 1(j≤2), . . . , 1(j≤I))0ξ= (µ11,∆α2, . . . , ∆αI,∆β2, . . . , ∆βI)0∈Rp. For the asymptotic theory in the over-dispersed Poisson model, it turns out to be useful to explicitly decouple the level and its deviations by decomposing as: ξ= (µ11,ξ(2)0)0and xij = (1, x(2)0 ij )0. We can then define the aggregate predictor τand the frequencies πij as: τ=∑ ij∈I exp(µij)and πij =exp(µij) τ=exp(x(2)0 ij ξ(2)) ∑ij∈I exp(x(2)0 ij ξ(2)). (2) Importantly, the frequencies πij are invariant to the level µ11 , the first component of ξ . Therefore, we can vary the aggregate predictor τ by varying µ11 without affecting the frequencies πij . The frequencies πij are, in turn, functions of ξ(2) alone. Further, we note that given ξ(2) , there is a one-to-one mapping between µ11 and τthrough τ=exp(µ11)∑ij∈I exp(x(2)0 ij ξ(2)). While this choice of identification scheme is useful for derivation of the theory in this paper, any scheme may be used in applications of the results. This is because, as Kuang et al. (2008) point out, the linear predictor µij is identified, unlike the individual effects. Since the main results of the paper rely on estimates of the linear predictors µij alone, they are unaffected by the choice of a particular identification scheme. Furthermore, the results in this paper are not limited to the chain-ladder predictor; we could, for example, include a calendar effect. Nielsen (2015) derives the form of the design vector for extended chain-ladder predictors in generalized trapezoids. The identification method is implemented in the R (R Core Team 2017) package apc (Nielsen 2015), as well as in the Python package of the same name (Harnau 2017). We note that the identification method can introduce arbitrariness into the forecast for models that require parameter extrapolation, such as the extended chain-ladder model with calendar effects. In the standard chain-ladder model, we can forecast claim reserves without parameter extrapolation; in a continuous setting, Lee et al. (2015) refer to this as in-sample forecasting. In contrast, in the extended chain-ladder model, we cannot estimate parameters for future calendar years from the Risks 2018,6, 70 7 of 37 run-off triangle. For this case, Kuang et al. (2008) and Nielsen and Nielsen (2014) explain how forecasts can be influenced by ad hoc constraints and lay out conditions for the identification method that make forecasts invariant to these arbitrary and untestable constraints. 3.3. Over-Dispersed Poisson Model We give the assumptions of the over-dispersed Poisson model and discuss its estimation by Poisson quasi-likelihood. We state the sampling scheme proposed by Harnau and Nielsen (2017) and the asymptotic distribution of the estimators. 3.3.1. Assumptions The first assumption imposes the over-dispersed Poisson structure on the moments. We can write it as: MODP :E(Yij) = exp(µij),var(Yij) E(Yij)=σ2. The second assumption is distributional and allows for the asymptotic theory later on. We assume that the independent aggregate claims Yij have a non-degenerate, non-negative and infinitely-divisible distribution with at least three moments. As noted by Harnau and Nielsen (2017), an appealing example for claim reserving of such a distribution is compound Poisson. The interpretation is that the aggregate incremental claims Yij can be written as Yij =∑Nij `=1X` for a Poisson number of claims Nij D =Poisson{exp(µij)} independent of the independent and identically distributed random claim amounts X`. 3.3.2. Estimation We estimate the over-dispersed Poisson model by Poisson quasi-likelihood. The appeal is that, as noted in Section 3.2, Poisson quasi-likelihood estimation replicates the chain-ladder technique. We explicitly distinguish between the model, subscripted with ODP , and its standard estimators, subor super-scripted with ql , to avoid confusion later on when we evaluate the estimators under the rival model. The fitted values for the linear predictors are given by: ˆ µql ij =x0 ij ˆ ξql where ˆ ξql =arg max ξ∈Rp∑ ij∈I{Yij(x0 ijξ)−exp(x0 ijξ)}. The fitted value for the aggregate predictor τis then given by: ˆ τql =∑ ij∈I exp(ˆ µql ij ) = ∑ ij∈I Yij, a result implied by the fact that the re-parametrization of the Poisson likelihood in terms of the mixed parameter (τ , ξ(2)) is linearly separable, so the parameters are variation independent; see, for example, Harnau and Nielsen (2017); Martínez Miranda et al. (2015) or, for a more formal treatment, Barndorff-Nielsen (1978, Theorem 8.4). This implies that the estimator for the aggregate predictor is unbiased for the aggregate mean. As an estimator for the over-dispersion σ2,Harnau and Nielsen (2017) use the Poisson deviance D scaled by the degrees of freedom. The deviance is the log likelihood ratio statistic against a model with as many parameters as observations, giving a perfect fit. The estimator is given by: ˆ σ2=D n−pwhere D=2∑ ij∈I Yij{log(Yij)−ˆ µql ij }. (3) Risks 2018,6, 70 8 of 37 3.3.3. Sampling Scheme For the asymptotic theory, we adopt the sampling scheme proposed by Harnau and Nielsen (2017). The idea is to grow the overall mean τ=∑ij∈I E(Yij) while holding the frequencies πij and thus ξ(2) fixed. We note that this also implies that µ11 is O{log(τ)} . In this sampling scheme, information accumulates in the estimated frequencies. In this sense, it is reminiscent of multinomial sampling as used, for example, by Martínez Miranda et al. (2015) in a Poisson model conditional on the data sum. Furthermore, we assume that τ increases in such a way that the skewness vanishes. Harnau and Nielsen (2017) remark that this is implicit for distributions such as Poisson, negative binomial and many compound Poisson distributions. Importantly, the sampling scheme holds the number of cells in the run-off triangle fixed. If we instead grew the dimension of the array, the number of parameters would also increase, thus making an asymptotic theory difficult. 3.3.4. Asymptotic Theory Based on the assumptions in Section 3.3.1 and the sampling scheme Section 3.3.3,Harnau and Nielsen (2017) derived the asymptotic distribution of the estimators. The theory hinges on Harnau and Nielsen (2017, Theorems 1, 2), which for our purposes can be formulated as: τ−1/2{Yij −exp(µij)}=τ1/2(Yij/τ−πij)D →N(0, σ2πij)and Yij τ P →πij. (4) An implication of the sampling scheme is that we cannot consistently estimate µ11 since the overall mean τ and thus the level µ11 grow. However, the remaining parameters ξ(2) are fixed and can be estimated in a consistent way. To ease notation, we define the design matrix X and the diagonal matrix of frequencies Πso: X={xij :(i,j)∈ I}0and Π=diag{πij :(i,j)∈ I}. (5) Harnau and Nielsen (2017, Lemma 1) derive the distribution of the estimator for the mean parameters in terms of the mixed parametrization (τ , ξ(2)0)0 . The advantage is that the two components of the mixed parameter are variation independent, so the covariance matrix featured in the asymptotic distribution is block-diagonal. This property turns out to be useful for example in the derivation of distribution forecasts. However, we opt to state the results in terms of the original parameterization by ξ to ease the analogy with the generalized log-normal model below. For our purposes, this does not complicate the theory. As a corollary to Harnau and Nielsen (2017, Lemma 1), we can then state the distribution of the quasi-likelihood estimator ˆ ξql as follows. All proofs are in Appendix A. Corollary 1. In the over-dispersed Poisson model Sections 3.3.1 and 3.3.3, √τ(ˆ ξql −ξ) = √τ((ˆ µ11 −µ11) ˆ ξ(2) ql −ξ(2))D →N{0, σ2(X0ΠX)−1}. Thus, even though the level µ11 →∞ , the difference between estimator and level (ˆ µ11 −µ11) vanishes in probability. We note that τX0ΠX corresponds to the Fisher information about ξ in a Poisson model. Further, Harnau and Nielsen (2017, Lemma 1) find that the asymptotic distribution of the deviance is proportional to a χ2: DD →σ2χ2 n−p. Thus, the estimator ˆ σ2has an asymptotic distribution, which is unbiased for σ2. Risks 2018,6, 70 15 of 37 and the variance to mean ratio are constant across cells. This manifests in the collapse of both R ODP and RGLN to a point mass at n. Further, Butler and Paolella (2008, Lemma 3) derive the support of R for a variety of cases depending on the properties of Aand B. Building on their work, we can prove the following result. Lemma 5. The distributions R GLN and R ODP have the same support. In non-degenerate cases, the support is (l,r)for 0<l<r<∞. The cumulative distribution functions and densities of ratios of quadratic forms admit saddle point approximations. We adapt the discussion in Butler and Paolella (2008) to our scenario in which ν=0; a setting that matches Lieberman (1994). We aim to approximate: P(R≤r) = Pe0Ae e0Be≤r=P(Xr≤0)where Xr=e0(A−rB)e First, we compute the eigenvalues of A−rB denoted λ1 , . . . , λn . We can write the cumulant generating function K(s) , the log of the moment generating function ϕ(s) = E[exp{s(A−rB)}] , of Xr and its `-th derivative as: K(s) = −1 2 n ∑ t=1 log(1−2sλt),K(`)(s) = ∂ ∂s` K(s) = (2`−2)!! n ∑ t=1λt 1−2sλt` where a!! =a(a− 2 )(a− 4 )··· is the double factorial with the usual definition that 0 !! = 1. The saddle point is the root: ˆ s:K(1)(ˆ s) = n ∑ t=1λt 1−2ˆ sλt=0. Except for the special case when all eigenvalues λt are zero, so K(1)(s) = 0, ˆ s is unique since K(1)(s) is strictly increasing. The former case occurs if and only if E(Xr) = 0, which is the case for r=trace(A)/trace(B). This case is dealt with separately. For the other cases, we compute: ˆ w=sgn(ˆ s)q−2K(ˆ s)and ˆ u=ˆ sqK(2)(ˆ s). Then, denoting by Φ( . ) and φ( . ) the standard normal cdf and density, respectively, the first order approximation to the cdf of R is: b P(R≤r) =    Φ(ˆ w) + φ(ˆ w)( ˆ w−1−ˆ u−1), if E(Xr)6=0 1 2+K(3)(0) 6√2πK(2)(0)3/2 , if E(Xr) = 0. This saddle point approximation is a special case of the more general form in Lugannani and Rice (1980). This is what Lieberman (1994) built on. Lugannani and Rice (1980) analysed the error behaviour for a sum of independent and identically distributed random variables and showed uniformity of the errors for a large sample. Butler and Paolella (2008) instead considered a fixed sample size and show uniformity of errors in the tail of the distribution. This seems appealing for our scenario, since we would expect the rejection region of the test to correspond to the tail of the distribution. 4.5. Power We show that the conjecture of a link between the power of the tests and variation in the means raised above in Section 4.1 is correct. To prove this, we consider a sequential asymptotic argument in which first, depending on the data generating process, τ becomes large or ω2 becomes small and then the means become “more dispersed” in a sense made precise below. Based on this argument, we can Risks 2018,6, 70 16 of 37 justify a one-sided test where the rejection region corresponds to the upper tail when the null model is generalized log-normal and to the lower tail when it is over-dispersed Poisson. The sequential asymptotics allows us to exclusively consider the impact of more dispersed means on R GLN and R ODP without worrying about the effect on the distribution of Rls , Rql , R∗ ls or R∗ ql . However, larger mean dispersion would be linked to changes in ξ(2) , a parameter that we keep fixed when deriving the asymptotic distribution of the test statistics in the first stage of the asymptotics. Therefore, we would expect the approximation quality achieved in the first stage to be affected by the second stage. The interpretation of the results is thus for a given first stage approximation quality, however large τor small ω2may be needed to achieve this. We model “more dispersed” means by increasing the variation in the frequencies πij and specifically by letting some frequencies go to zero. In this way, we do not make a statement about the means in absolute terms, but merely say that some cell means become large relative to others. For our analysis, we exclude cells for which estimation would yield a perfect fit; equivalently, we can impose that the frequencies do not exclusively vanish for perfectly-fitted cells. For example, in a chain-ladder model for the run-off triangle in Table 1, the corner cells ( 1, 10 ) and ( 10, 1 ) would be fit perfectly as they have their own parameters ∆β10 and ∆α10. To increase the variation in the frequencies πij , we decide on n−q cells of the run-off triangle for which we want the frequencies to vanish. We require that the remaining q cells with non-vanishing frequencies make up an array on which we can estimate a model with the same structure for the linear predictor µij as for the full data without obtaining a perfect fit. For example, for a chain-ladder model in which µij =αi+βj+δ , this would be the case for rectangular arrays with at least two columns and rows or for triangular arrays with at least three rows and columns. For the ease of notation, we sort rows and columns of the frequency matrix Π defined in (5) such that the cells with vanishing frequencies are in the bottom right block of the matrix. Then, for a q×q matrix Π1and an (n−q)×(n−q)matrix Π2, we define a new frequency matrix: Π(t)=s(t) Π10 0tΠ2!where s(t) = {trace(Π1) + t·trace(Π2)}−1(11) so s(t) takes care of the normalization such that the elements of Π(t) are still frequencies. The idea is to model the vanishing frequencies by letting t→ 0. Clearly, Π(1) corresponds to Π , whereas Π(0) has all frequencies in the the bottom right block equal to zero. We assume that Π16=q−1I , so that the limiting case does not correspond to a scenario without variation in the frequencies. Similarly, we sort rows and columns of the design matrix X to obtain a convenient partition. We sort the rows such that the q cells relating to non-vanishing frequencies are in the first q rows. Further, we sort the columns so the p1 , say, parameters relevant for these q cells are in the first p1 columns. Then, we can partition: X= X11 X12 X21 X22!(12) where X11 is q×p1 and X22 is (n−q)×(p−p1) . Column sorting ensures that X12 = 0. Imposing that there is no perfect fit for the q cells without vanishing frequencies implies that p1<q , so there are fewer parameters than cells. We are now interested in the properties of the large τ or small ω2 limiting distributions of the R statistics’ when some frequencies are small. For fixed t , Theorems 1and 2apply. Thus, for fixed frequencies Π(t) , the large τ and small ω2 distributions of the R statistics in an over-dispersed Poisson and generalized log-normal model, respectively, are: R(t) ODP =U0Π−1/2 (t)MΠ−1/2 (t)U U0M∗ (t)Uand R(t) GLN =U0MU U0Π1/2 (t)M∗ (t)Π1/2 (t)U Risks 2018,6, 70 17 of 37 where M∗ (t) is the weighted least squares orthogonal projection matrix I−X∗ (t)(X∗0 (t)X∗ (t))−1X∗0 (t) for X∗ (t)=Π1/2 (t)X. Thus, Π(t)enters not only directly, but also indirectly through M∗ (t). We study the tests’ power by looking at the limit of R (t) ODP and R (t) GLN as t→ 0. We reiterate that the sequential asymptotics neglect interactions between first stage asymptotics for large τ or small ω2 and the second stage small t asymptotics. A first intuition that neglects the potential influence of M∗ (t) may tell us that R (t) GLN should be well behaved while R (t) ODP blows up for small t . This turns out to be correct. Theorem 3. Let UD =N( 0, I) and let U1 contain the first q elements of U . Further define ˘ Π1=s( 0 )Π1 , ˘ X∗ 11 =˘ Π1/2 1X11 and ˘ M∗ 11 =I−˘ X∗ 11(˘ X∗0 11 ˘ X∗ 11)−1˘ X∗0 11. Then, as t →0, R(t) GLN a.s. →U0MU U0 1˘ Π1/2 1˘ M∗ 11 ˘ Π1/2 1U1 :=R(0) GLN while R(t) ODP a.s. →∞. Further, for α∈( 0, 1 ) , let q(t) GLN,α be the α -quantile of R (t) GLN and similarly for q(t) ODP,α . Then, R(t) ODP >q(t) GLN,αand R(t) GLN ≤q(t) ODP,αalmost surely as t →0. Theorem 3justifies one-sided tests and shows that the power of the tests under either data generating process goes to unity in the sequential asymptotic argument. Since the distribution of R ODP and R GLN coincides for equal means, the power of the tests to distinguish between the data generating processes comes entirely from the variation in means. As the mean variation becomes large, R ODP first order stochastic dominates R GLN . Thus, we can consider the lower tail of R ODP and the upper tail of R GLN as rejection regions. While still controlling the size of the test under the null, we gain power compared to two-sided tests as the mean variation increases. The denominator of R(0) GLN can be interpreted as “bad” weighted least squares residuals in a homoskedastic Gaussian model computed on just the subset of q cells with non-vanishing frequencies. For a brief intuition as to why only cells and parameters relating to X11 matter in the limit of the denominator, we consider weighted least squares estimation for Z=Xξ+Π−1/2 (t)e , taking Π(t) as given. We solve this by minimizing |Π1/2 (t)(Z−Xξ)|2 . For t> 0, the minimum is given by Z0Π1/2 (t)M∗ (t)Π1/2 (t)Z . When t= 0, the last n−q elements of Z and rows of X corresponding to the vanishing frequencies do not contribute to the norm. The same holds for the last p−p1 parameters in ξ that are then not identified. Thus, for t= 0, letting Z1 contain the first p1 elements of Z , the minimum of the norm equals Z0 1˘ Π1/2 1˘ M∗ 11 ˘ Π1/2 1Z1. 5. Simulations With the theoretical results for encompassing tests between over-dispersed Poisson and generalized log-normal models in place, we show that they perform well in a simulation study. First, we show that saddle point approximations to the limiting distributions R ODP and R GLN are very accurate. Second, we tackle an issue that disappears in the limit; namely, the choice between asymptotically identical estimators that generally differ in finite samples. We show that finite sample performance is indeed affected by this choice. However, we find that for some choices, finite sample and asymptotic distributions are very close. Third, we show that the tests have high power in finite samples and, considering the behaviour of the limiting distributions alone, that power increases quickly with the variation in means. For the simulations and empirical applications below, we use the Python packages quad_form_ratio (Harnau 2018b) and apc (Harnau 2017). The package was inspired by the R (R Core Team 2017) package apc (Nielsen 2015) with similar functionality. 5.1. Quality of Saddle Point Approximations We show that saddle point approximations work well compared to large Monte Carlo simulations. Risks 2018,6, 70 18 of 37 We consider three parameterizations. First, we let the design X correspond to that of a chain-ladder model for a ten-by-ten run-off triangle and set the frequency matrix Π to the least squares estimates b Πls =Π(ˆ ξ(2) ls ) of the Verrall et al. (2010) data in Table 1( VNJ ). Second, for the same design, we now set the frequency matrix to the least squares plug-in estimates based on a popular dataset by Taylor and Ashe (1983) (TA). We provide these data in the Appendix Ain Table A2. Third, we consider a design X for an extended chain-ladder model in an eleven-by-eleven run-off triangle and set Π to the least squares plug-in estimates of the Barnett and Zehnwirth (2000) data ( BZ ), also shown in the Appendix in Table A1. We remark that in the computations, we drop the corner cells of the triangles that would be fit perfectly in any case; this helps to avoid numerical issues without affecting the results. Given a data generating process R chosen from R ODP and R GLN , a design matrix X and a frequency matrix Π , we use a large Monte Carlo simulation as a benchmark for the saddle point approximation. First, we draw B= 10 7 realizations rb from R. For the Monte Carlo cdf b PMC( R ≤q) = B−1∑B b=1 1 (rb≤q) , we then find the quantiles qα , so b PMC( R ≤qα) = α for α= 0, 0.005, 0.01, . . . , 1. To compute the saddle point approximation b PSP( R ≤q) , we use the implementation of the procedure described in Section 4.4 in the package quad_form_ratio . Then, for each Monte Carlo quantile qα , we compute the difference b PSP( R ≤qα)−α . Taking the Monte Carlo cdf as the truth, we refer to this as the saddle point approximation error. Figure 1a shows the generalized log-normal saddle point approximation error b PSP( R GLN ≤ qα)−α plotted against α . One and two (pointwise) Monte Carlo standard errors pα(1−α)/B are shaded in blue and green, respectively. While the approximation errors for TA are generally not significantly different from zero, the same cannot be said for the other two sets of parameters. For the parameterizations VNJ and BZ , the errors start and end in zero and are negative in between. Despite statistically-significant differences, the approximation is very good with a maximum absolute approximation error of just over − 0.006. The errors in the tails are much smaller, as we might have expected given the results by Butler and Paolella (2008) discussed in Section 4.4. 0.00 0.25 0.50 0.75 1.00 α −0.006 −0.004 −0.002 0.000 BZ TA V NJ (a)b PSP(RGLN ≤qα)−α 0.00 0.25 0.50 0.75 1.00 α 0.000 0.001 0.002 0.003 0.004 0.005 BZ TA V NJ (b)b PSP(RODP ≤qα)−α Figure 1. Approximation error of the first order saddle point approximation to R GLN , shown in ( a ), and R ODP , displayed in ( b ). Monte Carlo simulation with 10 7 draws taken as the truth. One and two Monte Carlo standard errors shaded in blue and green, respectively. Figure 1b shows the plot for the approximation error to R ODP produced in the same way as Figure 1a. The approximation error is positive and generally significantly different from zero across parameterizations. Yet, the largest error is about 0.005 with smaller errors in the tails. We would argue that the saddle point approximation errors, while statistically significant, are negligible in applications. That is, using a saddle point approximation rather than a large Monte Carlo simulation is unlikely to affect the practitioner’s modelling decision. Risks 2018,6, 70 19 of 37 5.2. Finite Sample Approximations under the Null The asymptotic theory above left us without guidance on how to choose between test statistics R and estimators for the nuisance parameter Π that appears in the limiting distributions R. While the choice is irrelevant for large τ or small ω2 , we show that it matters in finite samples and that some combinations perform much better than others when it comes to approximation under the null hypothesis. In applications, we approximate the distribution of R by b R= R (b Π) . That is, defining the α quantile of b R as qb R α , we hope that P(R≤qb R α)≈α under the null hypothesis. To assess whether this is justified, we simulate the approximation quality across 16 asymptotically identical combinations of R -statistics and ratios of quadratic forms b R . We describe the simulation process in three stages. First, we explain how we set up the data generating processes for the generalized log-normal and over-dispersed Poisson model. Second, we lay out explicitly the combinations we consider. Third, we explain how we compute the approximation errors. As in Section 5.1, we point out that we drop the corner cells of the triangles in simulations. This aids numerical stability without affecting the results. For the generalized log-normal model, we simulate independent log-normal variables Yij , so log(Yij)D =N(x0 ijξ , ω2) . We consider three settings for the true parameters corresponding largely to the estimates from the same three datasets we used in Section 5.1, namely the Verrall et al. (2010) data ( VNJ ), Taylor and Ashe (1983) data ( TA ) and Barnett and Zehnwirth (2000) data ( BZ ). Specifically, we consider pairs (ξ , ω2) set to the estimated counterparts (ˆ ξls , ˆ ω2/s) for s= 1, 2. The estimates ˆ ω2 are 0.39 for VNJ , 0.12 for TA and 0.001 for BZ . Theory tells us that the approximation errors should decrease with ω2, thus as sincreases. For the over-dispersed Poisson model, we use a compound Poisson-gamma data generating process, largely following Harnau and Nielsen (2017) and Harnau (2018a). We simulate independent Yij =∑Nij `=1Xs where Nij D =Poisson{exp(x0 ijξ)} and Xs are independent Gamma distributed with scale σ2− 1 and shape (σ2− 1 )−1 . This satisfies the assumptions for the over-dispersed Poisson model in Sections 3.3.1 and 3.3.3. For the true parameters (τ , ξ(2) , σ2) , we consider three sets of estimates (sˆ τql , ˆ ξ(2) ls , ˆ σ2) from the same data as for the log-normal data generating process. We use least squares estimates ˆ ξ(2) ls so that the frequency matrix Π is identical within parameterization between the two data generating processes. The estimates for ˆ σ2 are 10,393 for VNJ , 52,862 for TA and 124 for BZ . Those for ˆ τql are 14, 633, 814 for VNJ , 34, 358, 090 for TA and 10, 221, 194 for BZ . Again, we consider s= 1, 2, but this time scaling the aggregate predictor. If this increases, so should the approximation quality. We recall that ξ(2) and τ pin down µ11 through the one-to-one mapping τ=exp(µ11)∑ij∈I exp(x(2)0 ij ξ(2)). Thus, multiplying τby scorresponds to adding log(s)to µ11. For a given data generating process, we independently draw B= 10 5 run-off triangles 4b={Yij,b:(i , j)∈ I} and compute a battery of statistics for each draw. First, we compute the four test statistics Rls , Rql , R∗ ls and R∗ ql as defined in (10) . Second, we compute the estimates for the frequency matrices Π based on least squares estimates, quasi-likelihood estimates and feasible weighted least squares estimates with least squares and with the quasi-likelihood first stage. This leads to four different approximations to the limiting distribution, which, dropping the subscript for the data generating process, we denote by: b Rls =R{Π(ˆ ξ(2) ls )},b Rql =R{Π(ˆ ξ(2) ql )},b R∗ ls =R{Π(ˆ ξ∗(2) ls )}, and b R∗ ql =R{Π(ˆ ξ∗(2) ql )}. Given a data generating process and a choice of test statistic and limiting distribution approximation (R , b R) , we approximate P(R≤qb R α) by Monte Carlo simulation. For each combination (R , b R) , we have B paired realizations; for example, Rb and the distribution b Rb are based on the triangle 4b . Denote the saddle point approximation to the cdf of b Rb as Gb(q) = b PSP(b Rb≤q) . Neglecting the saddle point approximation error, we then compute P(R≤qb R α) as b PMC(R≤qb R α) = Risks 2018,6, 70 20 of 37 B−1∑B b=1 1 {Gb(Rb)≤α} , exploiting that Gb(Rb)≤α whenever Rb≤G−1 b(α) = qb R α,b . We do this for α∈ A ={0.005, 0.01, . . . , 0.995}. To evaluate the performance, we consider three metrics: area under the curve of absolute errors (also roughly the mean absolute error), maximum absolute error and error at (one-sided) 5% critical values. We compute the area under the curve as AUC =∑199 `=1|b PMC(R≤qb R α`)−α`|∆α` where α`= 0.005 ·` , so ∆α`= 0.005; we can also roughly interpret this as the mean absolute error MAE = 200 / 199 ·AUC since α`= 200 −1 . The maximum absolute error is maxα∈A |b PMC(R≤qb R α)−α| . Finally, the error at 5% critical values is b PMC(R>qb R GLN,0.95)− 0.05 for the generalized log-normal and b PMC(R≤qb R ODP,0.05)−0.05 for the over-dispersed Poisson data generating process. Figure 2shows bar charts for the area under the curve for all 16 combinations of R and b R stacked across the three parameterizations for s= 1. The chart is ordered by the sum of errors across parameterizations and data generating processes within combination, increasing from top to bottom. The maximum absolute error summed over parameterizations is indicated by “ + ”. Since a bar chart for the maximum absolute errors is qualitatively very similar to the plot for the area under the curve, we do not discuss it separately and instead provide it as Figure A1 in the Appendix A. 0.00.10.20.3 (Rql,b R∗ ls) (Rql,b Rls) (Rql,b R∗ ql) (Rql,b Rql) (Rls,b R∗ ql) (Rls,b R∗ ls) (Rls,b Rls) (R∗ ql,b R∗ ql) (Rls,b Rql) (R∗ ql,b R∗ ls) (R∗ ql,b Rls) (R∗ ls,b Rql) (R∗ ql,b Rql) (R∗ ls,b R∗ ql) (R∗ ls,b Rls) (R∗ ls,b R∗ ls) Generalized Log-Normal 0.0 0.1 0.2 0.3 Over-dispersed Poisson V NJ TA BZ Figure 2. Bar chart of the area under the curve of absolute approximation errors (also roughly mean absolute error) for the considered combinations of R and b R . Ordered by the sum of errors within the combination across the data generating processes and parameterizations increasing from top to bottom. Sum of maximum absolute errors across parameterizations indicated by “ + ”. VNJ , TA and BZ are short for parameters set to their estimates from the Verrall et al. (2010) data in Table 1, the Taylor and Ashe (1983) data in Table A2 and the Barnett and Zehnwirth (2000) data in Table A1, respectively. Based on 105repetitions for each parametrization. s=1. Looking first at the sum over parameterizations and data generating processes within combinations, we see large differences in approximation quality both for the area under the curve of absolute errors and the maximum absolute error. The former varies from about 5 pp (percentage-points) for (R∗ ls , b R∗ ls) to close to 30 pp for (Rql , b R∗ ls) , the latter from 8 pp to 45 pp. It is notable that the four combinations involving Rql are congregated at the bottom of the pack. In contrast, the three best performing combinations all involve R∗ ls . These three top-performers have a substantial head start compared to their competition. While their AUC varies from 4.8 pp to 6.0 pp, there is a jump to 13.9 pp for fourth place. Similarly, the maximum absolute errors of the top three contenders lie between 7.5 pp and 9.3 pp, while those for fourth place add up to 20.5 pp. Risks 2018,6, 70 21 of 37 Considering next the contributions of the individual parameterizations to the area under the curve across data generating processes, the influence is by no means balanced. Instead, the average contribution over combinations of the VNJ , TA and BZ parameterizations is about 35%, 57% and 8%, respectively. This ordering is well aligned in magnitude and ordering with that of ω2 and σ2/τ , loosely interpretable as a measure for the expected approximation quality. Still, considering the contributions of the parameterizations within combinations, we see substantial heterogeneity. For example, the TA parameterization contributes much less to (R∗ ql,b R∗ ql)than VNJ, while the reverse is true for (Rls,b Rls). Finally, we see substantial variation between the two data generating processes. While the range of areas under the curve of absolute errors aggregated over parameterizations for the generalized log-normal is 0.7 pp to 10 pp, that for the over-dispersed Poisson is 3.2 pp to 19.4 pp. The best performer for the generalized log-normal is, perhaps unsurprisingly, (R∗ ls , b Rls) . Intuitively, since the data generating process is log-normal, the asymptotic results would be exact for this combination if we plugged the true parameters into the frequency matrices. Just shy of these, we plug in the least squares parameter estimates, which are maximum likelihood estimated. It is perhaps more surprising that using Rql is not generally a good idea for the over-dispersed Poisson data generating process even though the fact that these combinations take the bottom four slots is largely driven by the TA parametrization. Reassuringly, the top three performers across data generating processes also take the top three spots within data generating processes, albeit with a slightly changed ordering. Figure 3a shows box plots for the size error at 5% nominal size computed over the three parameterizations and two data generating processes within combinations (R , b R) for s= 1. Positive errors indicate an over-sized and negative errors an under-sized test. In the plots, medians are indicated by blue lines inside the boxes. The boxes show the interquartile range. Whiskers represent the full range. The ordering is increasing in the sum of the absolute errors at 5% critical values from top to bottom. 0.00 0.02 0.04 (Rql,b Rls) (Rql,b R∗ ls) (Rql,b R∗ ql) (R∗ ls,b Rql) (Rql,b Rql) (Rls,b Rls) (Rls,b R∗ ls) (Rls,b R∗ ql) (R∗ ql,b Rls) (R∗ ql,b Rql) (R∗ ql,b R∗ ql) (R∗ ql,b R∗ ls) (R∗ ls,b R∗ ql) (Rls,b Rql) (R∗ ls,b Rls) (R∗ ls,b R∗ ls) (a)s=1 0.00 0.02 0.04 (Rql,b Rls) (Rql,b R∗ ls) (Rql,b R∗ ql) (Rls,b Rls) (R∗ ql,b Rls) (Rls,b R∗ ls) (Rls,b R∗ ql) (R∗ ql,b R∗ ql) (R∗ ql,b R∗ ls) (Rql,b Rql) (R∗ ql,b Rql) (R∗ ls,b Rql) (R∗ ls,b R∗ ql) (Rls,b Rql) (R∗ ls,b Rls) (R∗ ls,b R∗ ls) (b)s=2 Figure 3. Box plots of size error at 5% critical values over parameterizations ( VNJ , TA and BZ ) and data generating processes (generalized log-normal and over-dispersed Poisson). Results for s= 1 shown in ( a ) and for s= 2 in ( b ). Medians indicated by blue lines inside the boxes. The boxes show the interquartile range. Whiskers represent the full range. Looking at the medians, we can see that these are close to zero, ranging from − 0.15 pp to 0.37 pp. However, there is substantial variation in the interquartile range, 0.1 pp for (Rls , b R∗ ql) to 0.7 pp for (R∗ ql , b Rls) and range, 0.4 pp for (R∗ ls , b R∗ ls) to 6.7 pp for (Rql , b Rls) ). The best and worst performers from the analysis for the area under the curve and maximum absolute errors are still found in the top and Risks 2018,6, 70 22 of 37 bottom positions. Particularly the performance of (R∗ ls , b R∗ ls) seems close to perfection with a range from −0.2 pp to 0.2 pp. Figure 3b is constructed in the same way as Figure 3a, but for s= 2, halving the variance for the generalized log-normal and doubling the aggregate predictor for the over-dispersed Poisson data generating process. Theory tells us that the approximation quality should improve, and this is indeed what we see. The medians move towards zero, now taking values between − 0.05 pp and 0.14 pp; the largest interquartile range is now 1.1 pp and the largest range 2.9 pp. Overall, the combination (R∗ ls , b R∗ ls) performs very well across the considered parameterizations and data generating processes. This is not to say that we could not marginally increase performance in certain cases, for example by picking (R∗ ls , b Rls) when the true data generating process is log-normal. However, even in this case in which we get the data generating process exactly right, not much seems to be gained in approximation quality where it matters most, namely in the tails relevant for testing. Thus, it seems reasonable to simply use (R∗ ls , b R∗ ls) regardless of the hypothesized model, at least for size control. 5.3. Power Having convinced ourselves that we can control size across a number of parameterizations, we show that the tests have good power. First, we consider how the power in finite sample approximations compares to power in the limiting distributions. Second, we investigate how power changes as the means become more dispersed based on the impact on the limiting distributions R GLN and R ODP alone, as discussed in Section 4.5. 5.3.1. Finite Sample Approximations Under the Alternative We show that combinations of R -statistics and approximate limiting distributions b R that do well for size control under the null hypothesis also do well when it comes to power at 5% critical values. The data generating processes are identical to those in Section 5.2 and so are the three considered parameterizations VNJ , BZ and TA . To avoid numerical issues, we again drop the perfectly-fitted corner cells of the triangles without affecting the results. To avoid confusion, we stress that we do not consider the impact of more dispersed means in this section. Thus, if we mention asymptotic results, we refer to large τ when the true data generating process is over-dispersed Poisson and for small ω2 when it is generalized log-normal, holding the frequency matrix Πfixed. For a given parametrization, we first find the asymptotic power. When the generalized log-normal model is the null hypothesis, we find the 5% critical values cGLN :P( R GLN >cGLN) = 0.05, using the true parameter values for Π . Then, we compute the power P( R ODP >cGLN) . Conversely, when the over-dispersed Poisson is the null model, we find cODP :P( R ODP ≤cODP) = 0.05 and compute the power P( R GLN ≤cODP) . Lacking closed-form solutions, we again use saddle point approximations, iteratively solving the equations for the critical values to a precision of 10−4. Next, we approximate the finite sample power of the top four combinations for size control in Section 5.2, (R∗ ls , b R∗ ls) , (R∗ ls , b Rls) , (Rls , b Rql) and (R∗ ls , b R∗ ql) , by the rejection frequencies under the alternative for s= 1. For example, say the generalized log-normal model is the null hypothesis, and we want to compute the power for the combination (R∗ ls , b R∗ ls) . Then, we first draw B= 10 5 triangles 4b from the over-dispersed Poisson data generating process. For each draw b , we find 5% critical values cb R GLN,ls,b:P(b R∗ GLN,ls,b>cb R GLN,ls,b) = 0.05. We compute these based on saddle point approximations, solving iteratively up to a precision of 10 −4 . Then, we approximate the power as B−1∑B b=1 1 {R∗ GLN,ls,b>cb R GLN,ls,b} . For the over-dispersed Poisson null hypothesis, we proceed equivalently, using the left tail instead. In this way, we approximate power for all three parameterizations and all four combinations. Before we proceed, we point out that we should be cautious to interpret power without taking into account the size error in finite samples. A test with larger than nominal size would generally have Risks 2018,6, 70 23 of 37 a power advantage purely due to the size error. One way to control for this is to consider size-adjusted power, which levels the playing field by using critical values not at the nominal, but at the true size. In our case, this would correspond to critical values from the true distribution of the test statistic R , rather than the approximated distribution b R . Therefore, the choice of b R would not play a role any more. To sidestep this issue, we take a different approach and compare how close the power of the finite sample approximations matches the asymptotic power. Table 2shows the asymptotic power and the gap between power in finite sample approximations and asymptotic power. Looking at the asymptotic power first, we can see little variation between data generating processes within parameterizations. The power is highest for the VNJ parameterization with 99%, followed by BZ with 95% and TA with 65%. This ordering aligns with that of the standard deviations of the frequencies πij under these parameterizations, which are given by 0.016, 0.012 and 0.009 for VNJ,BZ and TA, respectively. Table 2. Power in % at 5% critical values for large τ (over-dispersed Poisson DGP) and small ω2 (generalized log-normal DGP) along with the power gap in pp for the top four performers from Table 3. DGP is short for data generating process. Based on 105repetitions. s=1. P(RGLN ≤cODP)P(R≤cb R ODP)−P(RGLN ≤cODP) H0DGP Π(R∗ ls,b R∗ ls) (R∗ ls,b Rls) (Rls,b Rql) (R∗ ls,b R∗ ql) GLN ODP VNJ 99.02 0.25 0.14 0.18 0.22 BZ 94.61 −0.94 −0.96 −0.97 −0.94 TA 65.39 4.18 3.41 5.28 3.15 P(RODP >cGLN )P(R>cb R ODP)-P(RODP >cGLN ) (R∗ ls,b R∗ ls) (R∗ ls,b Rls) (Rls,b Rql) (R∗ ls,b R∗ ql) ODP GLN VNJ 99.23 −0.30 −0.25 −0.34 −0.20 BZ 94.67 −1.14 −1.15 −1.12 −1.12 TA 64.73 0.81 0.15 4.68 2.49 When considering the finite sample approximations, we see that their power is relatively close to the asymptotic power. For VNJ , absolute deviations range from 0.14 pp to 0.34 pp and for BZ from 0.94 pp and 1.15 pp. Compared to that, discrepancies for the TA parameterization are larger. The smallest discrepancy of 0.14 pp arises for (R∗ ls , b Rls) when the data generating process is generalized log-normal. As before, this is intuitive since it corresponds to plugging maximum likelihood estimated parameters ˆ ξ(2) into Π . With 5.28 pp, the largest discrepancy arises for (Rls , b Rql) for an over-dispersed Poisson data generating process. Mean absolute errors across parameterizations and data generating processes are rather close, ranging from 1.01 pp for (R∗ ls , b Rls) to 2.1 pp for (Rls , b Rql) . Our proposed favourite from above (R∗ ls , b R∗ ls) comes in second with 1.27 pp. We would argue that we can still justify the use of (R∗ ls,b R∗ ls)regardless of the data generating process. 5.3.2. Increasing Mean Dispersion in Limiting Distributions We consider the impact of more dispersed means on power based on the the test statistics’ limiting distributions R GLN and R ODP . We show that the power grows quickly as we move from identical means across cells to a scenario where a single frequency hits zero. For a given diagonal frequency matrix Πwith values πij, we define the linear combination: Π(t)=tΠ+ (1−t)n−1In. Thus, for t= 1, we recover Π , while for t= 0, we are in a setting where all cells have the same frequencies, so all means are identical. In the latter scenario, R GLN and R ODP collapse to a Risks 2018,6, 70 24 of 37 point-mass at n , as discussed in Section 4.4. We consider t ranging from just over zero to just under tmax :tmax minij∈I(πij) + ( 1 −tmax)n−1= 0. The significance of tmax is that Π(tmax) corresponds to the matrix where the smallest frequency is exactly zero. For each t , we approximate one-sided 5% critical values of R (t) GLN = R GLN(Π(t)) and R (t) ODP = RODP(Π(t))through: c(t) GLN :b PSP(R(t) GLN >c(t) GLN) = 0.05 and c(t) ODP :b PSP(R(t) ODP ≤c(t) ODP) = 0.05. We iteratively solve the equations up to a precision of 10 −4 . Theorem 3tells us that the critical values should grow for both models, but that c(t) GLN converges as t approaches tmax , while c(t) ODP goes to infinity. Then, for given t and critical values, we find the power when the null model is generalized log-normal P( R (t) ODP >c(t) GLN) and when the null model is over-dispersed Poisson P( R (t) GLN ≤c(t) ODP) . Again, we use saddle point approximation. Based on Theorem 3, we should see the power go to unity as tapproaches tmax. We consider the same parameterizations VNJ , TA and BZ of frequency matrices Π and design matrices X as above. The values for tmax are 1.083 for VNJ , 1.396 for TA and 1.103 for BZ . To avoid numerical issues, we again drop the perfectly-fitted corner cells from the triangles. In this case, while the power is not affected, the critical values are scaled down by the ratio of ˆ τls computed over the smaller array without corner cells to that computed over the full triangle. Since this is merely proportional, the results are not affected qualitatively. Figure 4a shows the power when the generalized log-normal model is the null hypothesis. For all considered parameterizations, this is close to 5% for t close to zero, increasing monotonically with t and approaching unity as tapproaches tmax, as expected. 0.0 0.5 1.0 t 0.0 0.2 0.4 0.6 0.8 1.0V NJ TA BZ (a)P(R(t) ODP >c(t) GLN) 0.0 0.5 1.0 t −0.02 −0.01 0.00 0.01 V NJ TA BZ (b)P(R(t) ODP >c(t) GLN)−P(R(t) GLN ≤c(t) ODP) Figure 4. Power as t increases from zero to tmax . Values for tmax are 1.083 for VNJ , 1.396 for TA and 1.103 for BZ . ( a ) shows power when the null model is generalized log-normal, and ( b ) shows the difference in power between the two models. For t= 1, where Π(t) corresponds to the least squares estimated frequencies from the data, the power matches what we found in Table 2. Figure 4b shows the difference in power between the two models plotted over t . For the three settings we consider, these curves have a similar shape and start and end at zero. Generally, the power is very comparable, with differences between − 2 pp and 1 pp, again matching our findings from Table 2for t=1. Risks 2018,6, 70 31 of 37 Johansen (1979, Theorems 7.1 and 7.3, Lemma 7.2), we showed that the mapping Y7→ ˆ ξql estimator is asymptotically equivalent to the linear mapping Y7→ (X0ΠX)−1X0Y. Meanwhile, the weighted least squares estimator maps Z7→ (X0ΠX)−1X0Π1/2Π1/2Z=ˆ ξ∗ . Thus, the only non-linear component of the mapping Y7→ ˆ ξ∗ is the transformation from Y to Z . However, while this mapping is non-linear in finite sample, for large τ it is equivalent to the linear map from Y to Π−1Y as seen in the proof of Lemma 1. Asymptotically, this conforms to sequentially applying the transformations Π−1 followed by (X0ΠX)−1X0Π1/2Π1/2 to Y . Taken together, the map reduces to (X0ΠX)−1X0 . Thus, both the Poisson quasi-likelihood and weighted least squares mapping asymptotically apply the same transformation to τ−1/2{Y−exp(µ)}D →N( 0, σ2Π) . Thus, √τ(ˆ ξ∗−ˆ ξql)P →0. The proof for τRSS∗−DP → 0 follows by the same argument. The main insight is that the asymptotic distribution of the Poisson deviance is asymptotically equivalent to that of the quadratic form τ−1{Y−exp(µ)}0(Π−1−X0(X0ΠX)−1X0){Y−exp(µ)} , as Harnau and Nielsen (2017, Proof of Lemma 1) show building on Johansen (1979, Theorems 7.7, 7.8). This is again asymptotically identical to the sequential mapping from Y to Z followed by the map from Z to the scaled residual sum of weighted least squares τRSS∗={√τΠ1/2(Z−µ)}0M∗{√τΠ1/2(Z−µ)}. To show that we can replace the weight matrix Π in the weighted least squares estimator by b Πls or b Πql we note that both matrices converge in probability to Π and then apply Slutsky’s theorem (Casella and Berger 2002, Theorem 5.5.17). Combining this argument with the proof of the equivalence of D and RSS∗ in the last paragraph, it also follows that we can replace the weights in RSS∗ without affecting the result. Finally, both ˆ τls/τP →1 and ˆ τql/τP →1 so we can replace τas well by Slutsky’s theorem. Appendix A.5. Proof of Theorem 1 Taking into account the results from Lemma 3, it follows that (τRSS∗)/DP → 1 and that the result still holds if we replace RSS∗ by RSS∗ ls or RSS∗ ql and τ by ˆ τls or ˆ τql . Thus, for example, Rls P →Rql so their difference vanishes and similarly for any other of the six total combinations. Both RSS =Z0MZ and RSS∗=Z∗0M∗Z∗=ZΠ1/2 M∗Π1/2Z are quadratic forms in the same random vector Z . It follows from the proofs of Lemma 2and Lemma 3that τRSS D → U0Π−1/2MΠ−1/2U and τRSS∗D →U0M∗U for the same UD =N( 0, In) . The distribution of RSS/RSS∗ follows by the continuous mapping theorem. Since τRSS∗−DP → 0 as in Lemma 3, τRSS∗/DP → 1 so that RSS/RSS∗−τRSS/DP →0 follows. We can replace τ by ˆ τql since Harnau and Nielsen (2017, Theorem 2) gives us that τ/ˆ τql P → 1. From Lemma 2, ˆ τls/τP → 1. Further, both Π(ˆ ξ(2) ls ) and Π(ˆ ξ(2) ql ) converge to Π in probability. Then, by Slutsky’s theorem, we can replace the true parameters with their estimates in τRSS/D and RSS/RSS∗ without affection the limiting distribution. Appendix A.6. Proof of Lemma 4 The asymptotic distribution of the weighted least squares estimators follows by the same argument as in Lemma 2, except now taking (ω2)−1/2(Z−µ)D →N( 0, In) , as shown by Kuang and Nielsen (2018, Theorem 3.3), as a starting point. The asymptotic equivalence of weighted least squares and Poisson quasi likelihood estimation follows from the same argument as in Lemma 3, except now (ω2)−1/2{Y−exp(µ)}D → N[ 0, diag{exp(µ)}] . The argument for replacing true parameters in the frequency matrix Π and the aggregate means τby estimates is identical to that in Lemma 3as well. Risks 2018,6, 70 32 of 37 Appendix A.7. Proof of Theorem 2 This follows by the same argument as the proof for Theorem 1above, except now combining the asymptotic distribution of the least squares estimator in the generalized log-normal model from Kuang and Nielsen (2018, Theorem 3.5) in (8) and Lemma 4. Appendix A.8. Proof of Lemma 5 First, we show that RGLN and RODP share a common support. If we recall that: RGLN(U) = U0MU U0Π1/2M∗Π1/2Uand RODP(U) = U0Π−1/2MΠ−1/2U U0M∗U, the main insight is that R GLN(Π−1/2U) = R ODP(U) . Formally, both R GLN :Rn7→ R and R ODP : Rn7→ R are random variables on {Rn , B(Rn) , P} where P is the measure associated with N( 0, In) . We now show that P{ R GLN ∈ S} = 1 implies that P{ R ODP ∈ S} = 1; the opposite direction follows. Notationally, P{ R GLN ∈ S} =P{u∈Rn: R GLN(u)∈ S} and, for some set A ⊆ R we denote by R −1 GLN(A) the pre-image {u∈Rn: R GLN(u)∈ A} . Now, R GLN is measurable since it is continuous almost everywhere, the exception being the measure zero set where the denominator is zero. Thus, P{RGLN ∈ S} =P{U∈R−1 GLN(S)}. Since the support of N( 0, In) is Rn , we must have that R −1 GLN(S) = Rn and hence S= R GLN(Rn) . Since Π is invertible, Π−1/2Rn=Rn so that R GLN(Rn) = R GLN(Π−1/2Rn) = R ODP(Rn) . Thus, RODP(Rn) = Sand so R−1 GLN(S) = R−1 ODP(S). Taken together, P{RGLN ∈ S} =P{U∈R−1 GLN(S)}=P{U∈R−1 ODP(S)}=P{RODP ∈ S} =1. Since this holds in both directions and for any such S, it holds for the support which is a special case of S. Now that we showed that R GLN and R ODP have identical support, we show that the support is a bounded compact set (l , r) . To do so, we specify the support of R GLN . The key insight is that the real symmetric matrices A=M and B=Π1/2M∗Π1/2 commute so AB =BA . Thus, they are simultaneously diagonalizable (Newcomb 1961). Beyond that, both matrices are of rank n−p . Thus, we can find an orthogonal matrix of Eigenvectors P0such that: PAP0= ΛA0 0 0!and PBP0= ΛB0 0 0! where ΛA and ΛB are diagonal n−p×n−p matrices of Eigenvalues. Since M is a projection matrix, ΛA=I . Then, making use of Butler and Paolella (2008, Lemma 3 Case 2(c)), the upper bound of the support of R GLN is given by the largest element of Λ−1 BΛA=Λ−1 B . Thus, it is finite. To find the lower bound, we consider the upper bound of − R GLN , thus swapping A for −A . By the same argument, the lower bound is the largest element of Λ−1 B times − 1, thus the smallest element of Λ−1 B . Since B is positive semi-definite and because ΛB contains only the non-zero eigenvalues, this must be larger then zero. Denoting the diagonal elements of ΛB sorted in descending magnitude by λB= (λB,(1), . . . λB,(n−p)), we can write the support as (l, r) = (λ−1 B,(1),λ−1 B,(n−p)). Appendix A.9. Proof of Theorem 3 We assume that rows and columns from the design X and the frequency matrix Π that relate to cells with perfect fit, such as the corners in a run-off triangle, have been removed. In this way, there is no need to keep track of the restriction that not only frequencies relating to such cells can go to zero. Risks 2018,6, 70 33 of 37 Further, we assume without loss of generality that rows and columns of X and Π have been sorted as described in Section 4.5. First, we establish a useful equivalence of the denominators of both R (t) GLN and R (t) ODP . We can write these as a quadratic form of the type V(t)0M∗ (t)V(t) where V(t)=Π1/2 (t)U for R (t) GLN and V(t)=U for R(t) ODP. Now, for any t>0, V(t)0M∗ (t)V(t)=|M∗ (t)V(t)|2=min ξ∈Rp|V(t)−X∗ (t)ξ|2. To understand the last equality, we note that minimizing argument for the last expression is the least squares estimator ˆ ξ= (X∗0 (t)X∗ (t))−1X∗0 (t)V(t) . The expression itself corresponds to the squared length of the least squares residuals. We compute these residuals as M∗ (t)V(t) and their squared length is given by |M∗ (t)V(t)|2. Now, we derive the small t limit of R (t) GLN . We denote the first q and last n−q elements of U by U1 and U2 , respectively, and similarly for ξ . Then, using the partition of Π(t) in (11) and X in (12) where X12 =0, we can write the denominator as: U0Π1/2 (t)M∗Π1/2 (t)U=s(t)min ξ∈Rp Π1/2 1(U1−X11ξ1) tΠ1/2 2{U2−(X21ξ1+X22ξ2)}! 2 where s(t) = {trace(Π1) + t·trace(Π2)}−1 . We note that the norm can only every decrease as t becomes smaller. This is because for t1<t2 , we could still choose the optimal ξ∗ under t2 and make the norm smaller by replacing t2with t1. For t=0, the right hand side simplifies to: min ξ1∈Rp1|˘ Π1/2 1(U1−X11ξ1)|2=U0 1˘ Π1/2 1˘ M∗ 11 ˘ Π1/2 1U1. Therefore, lim t→0U0Π1/2 (t)M∗ (t)Π1/2 (t)U=U0 1˘ Π1/2 1˘ M∗ 11 ˘ Π1/2 1U1. The numerator does not depend on t. Thus, as long as U6=0, which happens with probability one, lim t→0R(t) GLN =U0MU limt→0U0Π1/2 (t)M∗ (t)Π1/2 (t)U=U0MU U0 1˘ Π1/2 1˘ M∗ 11 ˘ Π1/2 1U1 =R(0) GLN. It follows that R(t) GLN a.s. →R(0) GLN. Next, we consider R(t) ODP. For the denominator, we have: U0M∗ (t)U=min ξ∈Rp|U−X∗ (t)ξ|2≤ |U|2 where the inequality follows since we can always set ξ= 0. Thus, the denominator is bounded from above. Looking at the numerator U0Π−1/2 (t)MΠ−1/2 (t)U , we see that the limit of Π−1/2 (t) does not exist. We look at this in more detail. Partition M so M11 is the q×q top right element and so on. Then we can write the numerator as: s(t)−1(U0 1Π−1/2 1M11Π−1/2 1U1+t−1U2Π−1/2 2M22Π−1/2 2U2+2t−1/2U1Π−1/2 1M12Π−1/2 2U2). The normalization s(t) converges to {trace(Π1)}−1 as t→ 0. The first term of the sum is non-negative and does not vary with t . The second term is non-negative and Op(t−1) while the third term is Op(t−1/2) with ambiguous sign. Further, dropping cells with perfect fit from X ensures Risks 2018,6, 70 34 of 37 M22 6= 0. Thus, since U26= 0 with probability one, the second term is positive with probability one so overall the numerator U0Π−1/2 (t)MΠ−1/2 (t)Ua.s. →∞ . Thus, since the denominator is bounded, R (t) ODP a.s. →∞ . Finally, R (t) ODP >q(t) GLN,α almost surely since for α∈( 0, 1 ) , the quantile q(t) GLN,α a.s →q(0) GLN,α<∞ . Conversely, q(t) ODP,α a.s →∞so R(t) GLN ≤q(t) ODP,αalmost surely. 0.00.10.20.3 (Rql,b R∗ ls) (Rql,b Rls) (Rql,b R∗ ql) (Rql,b Rql) (Rls,b R∗ ls) (Rls,b R∗ ql) (Rls,b Rql) (Rls,b Rls) (R∗ ql,b R∗ ql) (R∗ ls,b Rql) (R∗ ql,b R∗ ls) (R∗ ql,b Rls) (R∗ ql,b Rql) (R∗ ls,b R∗ ql) (R∗ ls,b R∗ ls) (R∗ ls,b Rls) Generalized Log-Normal 0.0 0.1 0.2 0.3 Over-dispersed Poisson V NJ TA BZ Figure A1. Bar chart of maximum absolute errors for the considered combinations of R and b R . Ordered by the sum of errors within combination across data generating processes and parameterizations increasing from top to bottom. Sum of maximum absolute errors across parameterizations indicated by “ + ”. VNJ , TA , and BZ is short for parameters set to their estimates from the Verrall et al. (2010), Taylor and Ashe (1983) and Barnett and Zehnwirth (2000) data, respectively. Based on 10 5 repetitions for each parametrization. s=1. Table A1. Insurance run-off triangle taken from Barnett and Zehnwirth (2000, Table 3.5) as used in the empirical application in Section 6.2 and the simulations in Section 5. i, j 1 2 3 4 5 6 7 8 9 10 11 1 153,638 188,412 134,534 87,456 60,348 42,404 31,238 21,252 16,622 14,440 12,200 2 178,536 226,412 158,894 104,686 71,448 47,990 35,576 24,818 22,662 18,000 - 3 210,172 259,168 188,388 123,074 83,380 56,086 38,496 33,768 27,400 - - 4 211,448 253,482 183,370 131,040 78,994 60,232 45,568 38,000 - - - 5 219,810 266,304 194,650 120,098 87,582 62,750 51,000 - - - - 6 205,654 252,746 177,506 129,522 96,786 82,400 - - - - - 7 197,716 255,408 194,648 142,328 105,600 - - - - - - 8 239,784 329,242 264,802 190,400 - - - - - - - 9 326,304 471,744 375,400 - - - - - - - - 10 420,778 590,400 - - - - - - - - - 11 496,200 - - - - - - - - - - Risks 2018,6, 70 35 of 37 Table A2. Insurance run-off triangle taken from Taylor and Ashe (1983) as used in the empirical application in Section 6.3 and the simulations in Section 5. i, j 1 2 3 4 5 6 7 8 9 10 1 357,848 766,940 610,542 482,940 527,326 574,398 146,342 139,950 227,229 67,948 2 352,118 884,021 933,894 1,183,289 445,745 320,996 527,804 266,172 425,046 - 3 290,507 1,001,799 926,219 1,016,654 750,816 146,923 495,992 280,405 - - 4 310,608 1,108,250 776,189 1,562,400 272,482 352,053 206,286 - - - 5 443,160 693,190 991,983 769,488 504,851 470,639 - - - - 6 396,132 937,085 847,498 805,037 705,960 - - - - - 7 440,832 847,631 1,131,398 1,063,269 - - - - - - 8 359,480 1,061,648 1,443,370 - - - - - - - 9 376,686 986,608 - - - - - - - - 10 344,014 - - - - - - - - - References Barndorff-Nielsen, Ole E. 1978. Information and Exponential Families. Chichester: Wiley. Barnett, Glen, and Ben Zehnwirth. 2000. Best estimates for reserves. Proceedings of the Casualty Actuarial Society LXXXVII: 245–321. Bartlett, Maurice S. 1937. Properties of Sufficiency and Statistical Tests. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 160: 268–82. [CrossRef] Beard, Robert E., Teivo Pentikäinen, and Erkki Pesonen. 1984. Risk Theory: The Stochastic Basis of Insurance, 3rd ed. London and New York: Chapman and Hall. Butler, Ronald W., and Marc S. Paolella. 2008. Uniform saddle point approximations for ratios of quadratic forms. Bernoulli 14: 140–54. [CrossRef] Casella, George, and Berger, Roger L. 2002. Statistical Inference, 2nd ed. Pacific Grove: Duxbury/Thomson Learning. Casualty Actuarial Society. 2011. Loss Reserving Data Pulled From Naic Schedule P. Available online: http: //www.casact.org/research/index.cfm?fa=loss_reserves_data (accessed on 8 July 2018). Chow, Gregory C. 1960. Tests of Equality Between Sets of Coefficients in Two Linear Regressions. Econometrica 28: 591–605. [CrossRef] Cox, David R. 1961. Tests of separate families of hypothesis. In Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability. Berkeley: University of California Press, pp. 105–23. Cox, David R. 1962. Further Results on Tests of Separate Families of Hypotheses. Journal of the Royal Statistical Society, Series B 24: 406–24. Durbin, James, and Geoffrey S. Watson. 1950. Testing for serial correlation in least squares regression. I. Biometrika 37: 409–28. [PubMed] Durbin, James, and Geoffrey S. Watson. 1951. Testing for Serial Correlation in Least Squares Regression. II. Biometrika 38: 159–77. [CrossRef] [PubMed] England, Peter, and Richard Verrall. 1999. Analytic and bootstrap estimates of prediction errors in claims reserving. Insurance: Mathematics and Economics 25: 281–93. [CrossRef] England, Peter D., and Richard J. Verrall. 2002. Stochastic Claims Reserving in General Insurance. British Actuarial Journal 8: 443–518. [CrossRef] England, Peter. 2002. Addendum to “Analytic and bootstrap estimates of prediction errors in claims reserving”. Insurance: Mathematics and Economics 31: 461–66. [CrossRef] Ermini, Luigi, and David F. Hendry. 2008. Log income vs. linear income: An application of the encompassing principle. Oxford Bulletin of Economics and Statistics 70: 807–27. [CrossRef] Harnau, Jonas, and Bent Nielsen. 2017. Over-dispersed age-period-cohort models. Journal of the American Statistical Association. [CrossRef] Harnau, Jonas 2017. apc. Available online: https://pypi.org/project/apc/ (accessed on 8 July 2018). Harnau, Jonas. 2018a. Misspecification Tests for Log-Normal and Over-Dispersed Poisson Chain-Ladder Models. Risks 6: 25. [CrossRef] Harnau, Jonas. 2018b. quad_form_ratio. Available online: https://pypi.org/project/quad-form-ratio/ (accessed on 8 July 2018). Risks 2018,6, 70 36 of 37 Hendry, David F., and Bent Nielsen. 2007. Econometric Modeling: A Likelihood Approach. Princeton: Princeton University Press. Hendry, David F., and Jean-Francois Richard. 1982. On the formulation of empirical models in dynamic econometrics. Journal of Econometrics 20: 3–33. [CrossRef] Johansen, Søren. 1979. Introduction to the Theory of Regular Exponential Families. Copenhagen: Institute of Mathematical Statistics, University of Copenhagen. Johnson, Norman Lloyd, Samuel Kotz, and N. Balakrishnan. 1995. Continuous Univariate Distributions Volume 1, 2nd ed. Chichester: Wiley. Kremer, Erhard. 1982. IBNR-claims and the two-way model of ANOVA. Scandinavian Actuarial Journal 1982: 47–55. [CrossRef] Kremer, Erhard. 1985. Einführung in die Versicherungsmathematik, 7th ed. Gottingen: Vandenhoeck & Ruprecht. Kuang, Di, and Bent Nielsen. 2018. Generalized Log-Normal Chain-Ladder. arXiv. [CrossRef] Kuang, Di, Bent Nielsen, and Jens P. Nielsen. 2008. Identification of the age-period-cohort model and the extended chain-ladder model. Biometrika 95: 979–86. [CrossRef] Kuang, Di, Bent Nielsen, and Jens Perch Nielsen. 2008. Forecasting with the age-period-cohort model and the extended chain-ladder model. Biometrika 95: 987–91. [CrossRef] Kuang, Di, Bent Nielsen, and Jens P. Nielsen. 2015. The geometric chain-ladder. Scandinavian Actuarial Journal 2015: 278–300. [CrossRef] Lancaster, Tony. 2000. The incidental parameter problem since 1948. Journal of Econometrics 95: 391–413. [CrossRef] Lee, Young K., Enno Mammen, Jens P. Nielsen, and Byeong U. Park. 2015. Asymptotics for in-sample density forecasting. Annals of Statistics 43: 620–51. [CrossRef] Lieberman, Offer. 1994. Saddle point approximation for the distribution of a ratio of quadratic forms in normal variables. Journal of the American Statistical Association 89: 924–28. [CrossRef] Lugannani, Robert, and Stephen Rice. 1980. Saddle point approximation for the distribution of the sum of independent random variables. Advances in Applied Probability 12: 475–90. Mack, Thomas. 1993. Distribution free calculation of the standard error of chain ladder reserve estimates. ASTIN Bulletin 23: 213–25. [CrossRef] Martínez Miranda, María Dolores, Bent Nielsen, and Jens Perch Nielsen. 2015. Inference and forecasting in the age-period-cohort model with unknown exposure with an application to mesothelioma mortality. Journal of the Royal Statistical Society: Series A (Statistics in Society) 178: 29–55. [CrossRef] Mizon, Grayham E., and Jean-Francois Richard. 1986. The Encompassing Principle and its Application to Testing Non-Nested Hypotheses. Econometrica 54: 657–78. [CrossRef] Newcomb, Robert W. 1961. On the simultaneous diagonalization of two semi-definite matrices. Quarterly of Applied Mathematics 19: 144–46. [CrossRef] Neyman, Jerzy, and Elizabeth L. Scott. 1948. Consistent estimates based on partially consistent observations. Econometrica 16: 1–32. [CrossRef] Nielsen, Bent, and Jens P. Nielsen. 2014. Identification and forecasting in mortality models. The Scientific World Journal 2014: 347043. [CrossRef] [PubMed] Nielsen, Bent. 2015. apc: An R Package for Age-Period-Cohort Analysis. The R Journal 7: 52–64. R Core Team. 2017. R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Taylor, Greg C., and Frank. R. Ashe. 1983. Second moments of estimates of outstanding claims. Journal of Econometrics 23: 37–61. [CrossRef] Thorin, Olof. 1977. On the infinite divisibility of the lognormal distribution. Scandinavian Actuarial Journal 1977: 121–48. [CrossRef] Verrall, Richard, Jens Perch Nielsen, and Anders Hedegaard Jessen. 2010. Prediction of RBNS and IBNR claims using claim amounts and claim counts. ASTIN Bulletin 40: 871–87. Verrall, Richard J. 1991. On the estimation of reserves from loglinear models. Insurance: Mathematics and Economics 10: 75–80. [CrossRef] Risks 2018,6, 70 37 of 37 Verrall, Richard J. 1994. Statistical methods for the chain ladder technique. Casualty Actuarial Society Forum 1: 393–446. Vuong, Quang H. 1989. Likelihood Ratio Tests for Model Selection and Non-Nested Hypotheses. Econometrica 57: 307–33. [CrossRef] c  2018 by the author. 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/).