A nonparametric visual test of mixed hazard models

Statistics & Operations Research Transactions SORT 37 (2) July-December 2013, 153-174 Statistics & Operations Research Transactions c Institut d’Estad´ ıstica de Catalunya [email protected] ISSN: 1696-2281 eISSN: 2013-8830 www.idescat.cat/sort/ A nonparametric visual test of mixed hazard models Jaap Spreeuw1, Jens Perch Nielsen2and Søren Fiig Jarner3 Abstract We consider mixed hazard models and introduce a new visual inspection technique capable of detecting the credibility of our model assumptions. Our technique is based on a transformed data approach, where the density of the transformed data should be close to the uniform distribution when our model assumptions are correct. To estimate the density on the transformed axis we take advantage of a recently defined local linear density estimator based on filtered data. We apply the method to national mortality data and show that it is capable of detecting signs of heterogeneity even in small data sets with substantial variability in observed death rates. MSC: 62F10, 62N01, 62N02, 62P05. Keywords: Mortality data, frailty models, visual inspection. 1. Introduction There is an increasing use of mortality models to answer a number of pension related questions. Mortality tables and their estimation have always been of importance while calculating appropriate prices of risk products depending on individuals’ survival. More recently, mortality models are being used in more complex models assessing the value of financial products incorporating survival in a variety of ways. Financial users of mortality models are therefore not only actuaries nowadays, but also investors looking for opportunities in survival bonds and other packages of survival risks. Different purposes of mortality models lead to different measures of quality. 1Faculty of Actuarial Science and Insurance, Cass Business School, City University London, 106 Bunhill Row, London, EC1Y 8TZ, UK. E-mail: [email protected] 2Faculty of Actuarial Science and Insurance, Cass Business School, City University London, 106 Bunhill Row, London, EC1Y 8TZ, UK. 3Danish Labour Market Supplementary Pension Fund, Kongens Vænge 8, 3400 Hillerød, Denmark. Received: October 2012 Accepted: July 2013 154 A nonparametric visual test of mixed hazard models In this paper we develop a visualization technique that seems useful for the individual assessment of the quality of a mortality model. One application we are thinking of is forecasting of mortalities that is a basic building block for the financial pricing of survival, but also a useful tool in asset liability management of pension portfolios. Typically, relatively simple parametric mortality models including calendar effects are used as starting point for mortality forecasts. The calendar effect is the explicit tool for the forecast and is often isolated and estimated through standard time series methodology. A perfect historical fit of the past is therefore not always what the mortality modeller is looking for. Often it is more important to have an overall good fit, without too systematic deviations giving reliable and meaningful forecasts. These latter objectives are not easy to generalize to some quantitative model that can be tested. Often simple mortality models are rejected, simply because mortality data often is nationwide and sufficiently abundant to inform relatively complex underlying parametric structures. Therefore, a test rejecting our simple model is often not what we want. We do know that our simple model is not accurate, we do not want an excessive fit, what we want is a good, intuitive and reliable forecast. When modelling mortality of a population, there is a variety of potentially suitable lifetime data models available. Potential models differ in levels of complexity and they try to capture different features of data. Specific parametric life tables combined with time series forecasts are omnipresent in the actuarial and demographic literature. The literature about parametric mortality projection has been developing rapidly in the last few years. Recent reviews of mainstream mortality forecasting models can be found in Cairns et al. (2009), Cairns et al. (2011), Dowd et al. (2010a,b) and Haberman and Renshaw (2011). Cairns et al. (2009) compare eight models on the basis of several desirable ex post qualitative properties (like model parsimony, transparency, possibility to generate sample paths, presence (or absence) of cohort effects and ability to achieve a nontrivial correlation structure) and quantitative criteria (consistency with historical data and robustness of parameter estimates). Six of these models are subject of subsequent investigation by Dowd et al. (2010a,b) and Cairns et al. (2011). These include the original Lee-Carter model (Lee and Carter, 1992), the basic age-period-cohort model by Renshaw and Haberman (2006), an alternative age-period-cohort model by Currie (2006), the original Cairns-Blake-Dowd model as launched in Cairns et al. (2006), and two extensions thereof. The six models are the subject of formal goodness-of-fit tests in Dowd et al. (2010a) and backtesting in Dowd et al. (2010b). Cairns et al. (2011) judges these models on the basis of ex ante qualitative aspects like biological reasonableness, plausibility of forecast levels of uncertainty in projections at several ages, and robustness of forecasts. In all these papers, the mortality data applied was confined to those of individuals aged 60 or above. Haberman and Renshaw (2011), concentrating on the key factors of life expectancy and annuity values, first conduct a detailed comparison of the several models at pensioner ages. Apart from the models in the above papers, they also consider special cases of the Renshaw and Haberman (2006) model in their study. Later on, they extend the age range and involve the model by Plat (2009) and several variants thereof. Jaap Spreeuw, Jens Perch Nielsen and Søren Fiig Jarner 155 The stability of the forecast depends crucially on the choice of the parametric form. Generally, a complex model with many parameters is not a good choice even though such models might be selected from classical mathematical statistical model selection designed for in-sample prediction. Models with many parameters generally fit data better than models with fewer parameters. On the other hand, a large number of parameters are harder to forecast than fewer parameters. Forecasting uncertainty increases dramatically with the number of parameters. Thus, to obtain reliable forecasts we want models which describe the key features of data with as few parameters as possible. The purpose of this paper is to introduce a visual diagnostic tool which can be used to guide us when choosing a parametric model. A good parametric model is a simple model without obvious systematic errors. That model could be chosen by the well informed statistician working with the particular mortality forecast application in mind. Our visual diagnostic tool will be just one helpful tool in the overall mathematical statistical toolbox. Our method is inspired from recent developments in extreme value estimation, where transformations of data give visual information on the quality of the distributional fit in the tail. This recent methodology has found its way into insurance pricing and also the related field of operational risk. For a comprehensive overview of this new transformation methodology in the latter context, see Bolanc´ e et al. (2012a). The transformation based method can for example compare the performance of several candidate models for a data set at hand. Assume we were told by an oracle what the exact true distribution is, then we would transform our data using this oracle information such that our transformed data would exactly originate from a uniform distribution. Now we do not have access to any oracles. However, if we take some estimated parametrically fitted survival distribution as defining our transformation, then any detectable deviance on the transformed scale from the uniform distribution implies deviances of the parametric distribution used in the transformation step from the underlying true distribution. Our methodology uses a nonparametric smooth kernel estimator on the transformed scale. One difficulty we meet here is that our data is classical survival data that is not independent identically distributed. We therefore use a recent local linear kernel density estimator – specifically the one of Nielsen et al. (2009) – that is adjusted for the truncation and censoring pattern we meet in our data. Comparison between different underlying suggested parametric models are carried out by first estimating these parametric models and then to investigate through visual inspection, whether the density of the transformed data indeed looks uniform. If the underlying parametric model under investigation would be true, the estimated density should be close to one over the unit interval. Therefore different underlying parametric models can be visualized and compared on the transformed scale. In principle, the densities could also be estimated and compared on the original scale. However, there are several visual and estimational advantages to working on the transformed scale. One of these is that our method makes maximal use of sparse and volatile data and is thus particularly well suited to explore how potential models describe the mortality at ad- 156 A nonparametric visual test of mixed hazard models vanced ages where exposure is invariably limited. We test our method using data from nations of different size: USA, United Kingdom, Denmark and Iceland. Although the main focus of our paper is to model human mortality, it is worthwhile mentioning that our methodology is applicable to any probability density model, whether it concerns human survival or not. 1.1. Mixed hazard models Frailty theory offers a possible explanation to the presence of an old-age mortality plateau. According to this theory populations are heterogeneous with some people being more frail, i.e. having a higher hazard rate, than other people. Since persons with high hazard rates tend to die sooner than persons with low hazard rates old age groups will be dominated by low frailty persons and this effect reduces the rate of increase at the population level. Frailty models were introduced in the demographic literature by Vaupel et al. (1979). In a multiplicative frailty model, an individual’s hazard rate consists of two parts, namely a certain standard intensity and a certain nonnegative random variable, the frailty, acting multiplicatively on the standard intensity. A Gompertz or Makeham specification is usually taken for the standard intensity, although sometimes a Weibull model can be seen. Frailty is usually assumed to follow a Gamma distribution, which is known to be mathematically very tractable. A few publications about frailty modelling appeared in the actuarial literature. Wang and Brown (1998) use the Gompertz-Gamma or Perks model to graduate mortality improvement factors in a Society of Actuaries’ Life Table. Butt and Haberman (2004) employ Generalized Linear Models to graduate mortality of insured lives. They consider three mixture models, namely i) Perks; ii) modified Perks, and iii) Gompertz-Inverse Gaussian. The authors conclude that the Perks model fits the data best. An overview of heterogeneity models in life insurance is given in Olivieri (2006), while Jones (1998) develops a multiple state model to measure the impact of frailty on the propensity to lapse a policy. Finally, Li et al. (2009) extend the Lee-Carter model by allowing for unobserved heterogeneity within a cell, determined by age and time. In this paper we illustrate our methodology in the one dimensional case. Most forecasting models operate with a multiplicative relationship between age effect and time effect. To visualize the fit of the age effect, one would then have to divide out the estimated time effect and vice versa to visualize the time effect only. We are happy to say that our paper – diffused in preliminary versions – already has inspired a number of other works in mathematical and computational statistics. It has for example been cited in the three recent papers G´ amiz-P´ erez et al. (2013a,b,c). Jaap Spreeuw, Jens Perch Nielsen and Søren Fiig Jarner 157 1.2. Outline The set-up of this paper is as follows. In Section 2 we present the visual inspection technique in detail. Both the continuous-time framework with transformed counting processes and the implementation with discrete data is discussed. Section 3 discusses frailty models in general and introduces the class of models we will be using. Section 4 presents the numerical application. For four countries varying significantly in size (United States, United Kingdom, Denmark and Iceland), one data set per country (female period 2006 from the Human Mortality Database) and three different frailty specifications, namely Gamma, Inverse Gaussian, and degenerate (no frailty), we show the estimates as well as the visual inspection technique. In particular, we give a thorough analysis of the mortality at advanced ages that can be extracted from the continuous graphs. Section 5 sets out a conclusion. 2. Visual inspection technique 2.1. Sampling scheme of the survival data Consider a data set with mortality statistics of nlives. Let for each of these nindividuals Yibe an exposure process with value one when the i’th individual is alive and under observation and let Nibe a counting process taking the value one if the i’th individual has died while under observation. Both Yiand Niare functions of the age x. Formally, we assume that Niis a one-dimensional counting process with respect to an increasing right continuous complete filtration Fx,x∈R+,i.e. one that obeys les conditions habituelles, see Andersen et al. (1993, p. 60). We model the intensity as λc i(x) = µθ(x)Yi(x), where θbelongs to the parameter space Θof the parameters determining the exact mortality and frailty. The estimator b θof θis derived from minimizing the log likelihood of Borgan (1984): l(θ) = n ∑ i=1Zlog{µθ(x)}dNi(x)− n ∑ i=1Zµθ(x)Yi(x)dx, that is maximized over the parameter space Θ. 2.2. Visual inspection by transformations Assume that some oracle has given us the true underlying c.d.f. Fθ. Then consider the transformed counting processes Ni=Ni◦F−1 θdefined on [0,1].If our oracle really had 158 A nonparametric visual test of mixed hazard models told us the truth, then Niwould have stochastic intensity λi(y) = α(y)Yi(y), where Yi(y) = YiF−1 θ(y)with α(y) = 1/(1−y)corresponding to the hazard of the uniform distribution with density f(y) = α(y)expZy 0−α(s)ds=1, for y∈[0,1]. Another more statistical term for oracle information is prior information. It is that type of information that is external to the data set at hand. In our application below our prior information will always be some parametric specification of the model and our oracle candidate for the true c.d.f will be Fb θ,where b θis the estimated parameter in the specified parametric model. If Fb θreally is a good description of the true c.d.f. F, then our data should be uniformly distributed after a transformation by Fb θ. To be able to inspect the credibility of our oracle information or prior information or parametric assumptions, we estimate the density fbased on the filtered survival data N1,Y1,...,Nn,Ynon [0,1]and see whether it looks flat. This density estimator should have good boundary correction because it is defined on the transformed axis [0,1].We suggest to use the natural weighted local linear density estimator of Nielsen et al. (2009): b f(y) = n ∑ i=1ZKy,b(y−s)Yi(s)b S(s)dNi(s), where Ky,b(y−s) = a2(y)−a1(y)(y−s) a0(y)a2(y)−{a1(y)}2Kb(y−s), with Kb(y−s) = 1 bK(y−s b),(1) and aj(y) = n ∑ i=1ZKb(y−s)(y−s)jYi(s)ds, Jaap Spreeuw, Jens Perch Nielsen and Søren Fiig Jarner 159 and b S(s) = ∏ t≤sn1−db Λ(t)o, being the Kaplan-Meier estimate of the survival function, with b Λ(s) = n ∑ i=1 s Z0nY(n)(t)o−1 dNi(t), where Y(n)(t) = ∑n i=1Yi(s). 2.3. Implementing with discrete data In most real life applications we only have discretized versions of the stochastic processes Yiand Niavailable. First we need to define the relevant discretized time points H1,...,HKand the corresponding differences hk=Hk−Hk−1for k∈{1,...,K}, with H0=0. We define HK=inft;Fb θ(t) = 1for any plausible survival function Fθ. Discretized data are often defined as occurrences and exposures. Let respectively Ok= n ∑ i=1ZHk Hk−1 dNi(x) and Ek= n ∑ i=1ZHk Hk−1 Yi(x)dx. Now assume that we only observe these discrete occurrences – the Ok’s – and exposures – the Ek’s. Then a natural approximation of the log likelihood function l(θ)above to our discrete observations would be ld(θ) = ∑ k{logµθ(H∗ k)}Ok−∑ k µθ(H∗ k)Ek, where H∗ k= (Hk−1+Hk)/2 . Now consider discretized time points on the axis transformed by Fb θ.Let Hk= F∗ b θ(Hk).hk=Hk−Hk−1and H∗ k=Hk−1+Hk/2 for k∈{1,...,K}.Note that HK=1. Also note that often the discrete time points are equidistant before the time transformation but not thereafter. 160 A nonparametric visual test of mixed hazard models On the transformed axis with time, the series H∗ 1,...,H∗ Kis transformed into H∗ 1,...,H∗ K. We will have occurrences Ok=Ok and exposures Ek=Ek∗hk/hk. Assume that we were given the true c.d.f. with very large risk exposures Ek. Then on the original axis Ok∼µb θ(H∗ k)Ekhkwhile on the transformed axis Ok=Ok∼αb θH∗ kEkhk. If the model were the true one, the hazard rates OkEkon the transformed axis would be equal to 11−H∗ k, and hence the density functions would be constant at 1. The local linear density estimator on the transformed axis will in the discrete case be defined as b fd(y) = ∑ k Kd,y,b(y−H∗ k)b St d(H∗ k)Ok,(2) where Kd,y,b(y−s) = a2,d(y)−a1,d(y)(y−s) a0,d(y)a2,d(y)−{a1,d(y)}2Kb(y−s), aj,d(y) = K ∑ k=1 Kb(y−H∗ k)(y−H∗ k)jEk and b St d(H∗ k) = 0.5nb St d(Hk−1)+ b St d(Hk)o=0.5"exp(− k−1 ∑ i=1 hi Oi Ei)+exp(− k ∑ i=1 hi Oi Ei)#. The choice of the bandwidth bdepends on the availability of data. Large countries have a large risk exposure; then most of the deviation between the density estimate and 1 can be attributed to model uncertainty. In such cases, no or hardly any smoothing is required and bcan be small. For not so densely populated countries with small risk exposure, on the other hand, proper smoothing – with a larger bandwidth – is needed to compensate for parameter uncertainty. Jaap Spreeuw, Jens Perch Nielsen and Søren Fiig Jarner 161 3. Mixed hazard models In an individual frailty model the individual effect for a life’s mortality acts multiplicatively. Assume that a cohort consists of nindividuals. Then for the ith person of the cohort, the individual effect is represented by the random variable Ziand the conditional force of mortality at age x, given Zi=zi, is given by µ(x,zi) = ziµ(x),i∈{1,...,n}, with µ(x)denoting the standard force of mortality at age x– which is the force of mortality of a life with frailty level 1 – and all Ziindependent and identically distributed, with a mean equal to 1. In this paper we will assume that the individual hazard is of the form µ(x) = exp(a0+a1x+a2x2).(3) In the notation of Forfar et al. (1988) this model is labelled GM(0,3). Note that the special case a2=0 leads to the Gompertz model (GM(0,2)). The structure in (3) forms the basis for national and international mortality modelling in Jarner and Kryger (2011). We have dµ(x)/dx =µ(x)(a1+2a2x). It is reasonable to assume that mortality is increasing as a function of age. This would imply a1≥0 and a2≥0. Nonnegative estimates of a1and a2are also obtained in Jarner and Kryger (2011). The relative change of mortality as a function of age x– defined in Horiuchi and Coale (1990) as k(x) = dlnµ(x)dx – is a linear function of age: k(x) = a1+2a2x. The cohort mortality at age xis given as µθ(x) = E[Z|x]·µ(x), where E[Z|x] denotes the mean frailty of lives surviving to age x. Let LZdenote the Laplace transform of frailty at birth, i.e. LZ(s) = E[exp(−sZ)]. It can then be shown, see e.g. Hougaard (1984), that E[Z|x] = −L′ Z[M(x)] LZ[M(x)] with M(x) = x Z0 µ(s)ds. Hence the cohort mortality can be easily calculated for all frailty specifications with known Laplace transform. In the literature, the Gamma distribution has been by far the most popular specification in the frailty model. This is partly due to its mathematical tractability. Abbring and Van den Berg (2007) show that, under mild conditions 168 A nonparametric visual test of mixed hazard models 0.0 0.2 0.4 0.6 0.8 1.0 0.96 0.98 1.00 1.02 1.04 Transformedaxis Figure 6: Denmark: Local linear density estimator as in (2), with b =1/6, on the transformed scale: Gamma frailty (solid) compared with Inverse Gaussian (dotted) and no frailty (dashed). 0.0 0.2 0.4 0.6 0.8 1.0 0.80 0.85 0.90 0.95 1.00 1.05 1.10 Transformedaxis Figure 7: Iceland: Local linear density estimator as in (2), with b =1/2, on the transformed scale: Gamma frailty (solid) compared with Inverse Gaussian (dotted) and no frailty (dashed). Jaap Spreeuw, Jens Perch Nielsen and Søren Fiig Jarner 169 Therefore, there is no reason to fear that we have smoothed too much and that should be the reason for the small deviance. This USA study gives us some confidence that the Gamma frailty survival model is working well also for smaller data sets, where bigger fluctuations are to be expected. For the United Kingdom the Gamma frailty survival model also fits relatively well, but now with deviances up to five percent. Surprisingly the Danish Gamma frailty survival density has very small deviances with the biggest being less than three percent. Iceland is another case, deviances up to 20% are found and the two frailty models do not seem to improve the fit compared to having no frailty at all. Overall the conclusion from the graphs is that the Gamma frailty makes the best fit, the Inverse Gaussian less so, but with both frailty models being superior to having no frailty at all. If we take a closer look at the tail of the three fitted Danish survival models at the transformed scale, we can get some further insight into the question posed in the introduction. It is indeed very clear that the flattening out of the Gamma frailty density in the tail helps the fit. The Gamma frailty version is much closer to one around the tail with about half the deviance from one compared to the no-frailty density version. In general, the performance of Inverse Gaussian is somewhat between that of Gamma and no frailty. For Iceland, the curves are almost identical, due to the small estimate of σ2 and very similar estimates of the other parameters. In other words, for Iceland, the cases of Gamma frailty, Inverse Gaussian frailty and no frailty are nearly the same. For US and UK the Gamma specification clearly provides the best description of data of the candidates considered. Both the Inverse Gaussian and no frailty alternative deviate substantially more in the right tail than Gamma frailty. These two specifications both overestimate old age mortality substantially, while the Gamma frailty seems to capture the old age mortality plateau evident in data. Moreover, the right tail deviations of Gamma frailty is of the same magnitude as deviations for younger age segments, while the right tail deviations for Inverse Gaussian and no frailty seems to diverge. While the picture is less clear for Denmark, the frailty densities also here improve the description of old age mortality. It also seems that without frailty the deviation diverges in the right tail, but the magnitude of deviation is much smaller than for US and UK. In contrast to US and UK, the Gamma and Inverse Gaussian essentially perform equally well. Thus we conclude that there is enough information in data to indicate the presence of heterogeneity, but not enough information to distinguish between the different kinds of heterogeneity. Lastly, Iceland has so little exposure and so much uncertainty in data that even with the method derived in this paper we cannot distinguish between the models. A critical part of the study concerns the performance of the estimator for advanced ages. To this end, for each country we calculate the second largest and largest points of intersection of the estimator with the horizontal line (i.e. the two largest roots of the equation b fd(y) = 1) and translate this back into the corresponding ages. For comparability, we have left out in this investigation the late spike of no frailty and the Inverse Gaussian in the USA case. The results are given in Table 2 below. It is quite clear that Gamma frailty densities in all case are having the last crossing point. This indicates 170 A nonparametric visual test of mixed hazard models that the Gamma frailty density provides the best description of old age mortality among the considered models. In the USA case, we get almost to the age of 100 before our transformed density drops below one. The lowest last crossing for the gamma frailty is still quite high, namely 93 years. Above this last crossing point on the transformed scale, all the fitted parametric models seem to have too low densities. Thus, above the last crossing point our parametric models are overstating the possibility of dying. In other words, above the last crossing point all models seem to be on the safe side. In particular the Gamma frailty seems well behaved for annuity purposes. The density of very old are a bit too high, but rarely more than two percent, and these two percent are on the safe side when calculating for example annuities. Most of the extra old age mass is taken from the interval between the next last crossing point and the last crossing point, where the underlying parametric densities are overestimated in all cases. Therefore, while none of the densities are making a perfect fit, the Gamma frailty density is very close and with good properties for the annuity forecaster. It is on the safe side for the very old ages, with an overall annuity that seems to be close to the truth, overestimating the density in the very old ages, but compensating for that overestimation in an interval leading up to those old ages. Without frailty the deviations for old ages are substantially larger than with (gamma) frailty. The transformed density is below 1 which indicates that the probability of dying old is overestimated. At first glance this appears to be at odds with the fact that without frailty the old-age hazard is overestimated, cf. Figure 1. The explanation is that while the old-age hazard is overestimated the hazard is underestimated in the age groups below and therefore too many attain the (high) age of 90, say, after which they die too quickly. The model without frailty is on the safe when setting aside reserves for annuities for 40 year-olds, but if we were to use the same model for older age groups it would only be conservative up to a certain point. This clearly is not a desirable feature, and it illustrates the point that overstating the probability of dying old for one cohort is not necessarily a conservative assumption for other cohorts. Notice that it would be hard to get this kind of detailed information from testing the underlying densities or even from graphical visualization techniques on the original scale. Therefore, our simple transformation technique has enabled us to comfort the statistician forecasting mortality models based on simple underlying parametric survival distributions. Table 2: Second largest and largest crossing point of density estimator with horizontal line at 1. Country Gamma Inverse Gaussian No frailty Second largest Largest Second largest Largest Second largest Largest crossing crossing crossing crossing crossing crossing point point point point point point US 92.37 98.68 92.14 95.21 92.14 95.21 UK 89.51 96.85 89.15 95.34 87.61 92.10 Denmark 85.53 94.77 85.47 94.67 84.91 93.63 Iceland 84.38 93.02 84.34 93.01 84.34 93.01 Jaap Spreeuw, Jens Perch Nielsen and Søren Fiig Jarner 171 5. Conclusions We have developed a new visual inspection technique of survival models. It generalizes developments of transformation techniques of i.i.d. data, see for example Bolanc´ e et al. (2008, 2012a, 2012b, 2013). The method seems useful in many versions of follow-up studies, see for example Guill´ en et al. (2012) and Pinquet et al. (2011). We imagine it to be useful when the applied statistician wants the data to guide his intuition. The working methodology could be through running the knowledge loop cycle: Data→Visualization→New Assumption a number of times until the final assumptions seem intuitively reasonable and well behaved also according to more standard statistical techniques. All the mortality projection models discussed in the Introduction involve both an age and a time dimension. As mentioned in the Introduction one can use our onedimensional visualization technique for the age effect after having adjusted for the time effect and vice versa when visualizing the time effect. A full multidimensional version of our methodology is also possible. One could use multidimensional density estimation of filtered data to introduce a similar visual inspection technique to assessing the quality of mortality depending on both age and time. See for example Buch-Kromann and Nielsen (2012) for a recent multivariate density estimator that could be used in our visual diagnostic step after having transformed our data with our favourite forecasting mortality model. Transformations and visual fitting as developed in this paper would also seem relevant in other areas of actuarial science as, for example, reserving, see the recent papers Mart´ ınez-Miranda et al. (2012) and Kuang et al. (2011). Acknowledgement This project was funded by a research grant from The Actuarial Foundation and the Society of Actuaries. References Abbring, J. H. and Van den Berg, G. J. (2007). The unobserved heterogeneity distribution in duration analysis. Biometrika, 94 (1), 87–99. Andersen, P. K., Borgan, O., Gill, R. D. and Keiding, N. (1993). Statistical Models Based on Counting Processes. Springer-Verlag, New York. Bolanc´ e, C., Guillen, M. and Nielsen, J. P. (2008). Inverse beta transformation in kernel density estimation. Statistics and Probability Letters, 78, 1757–1764. Bolanc´ e, C., Guill´ en, M., Nielsen, J. P. and Gustafsson, J. (2012a). Quantitative Operational Risk Models. Chapman and Hall/CRC Finance Series, New York. Bolanc´ e, C., Ayuso, M. and Guill´ en, M. (2012b). A nonparametric approach to analysing operational risk with an application to insurance fraud. The Journal of Operational Risk, 7 (1), 1–16. 172 A nonparametric visual test of mixed hazard models Bolanc´ e, C., Guill´ en, M., Gustafsson, J. and Nielsen, J. P. (2013). Adding prior knowledge to quantitative operational risk models. The Journal of Operational Risk, 8 (1), 17–32. Borgan, O. (1984). Maximum likelihood estimation in parametric counting process models, with applications to censored failure time data. Scandinavian Journal of Statistics, 11, 1–16. Buch-Kromann T. and Nielsen, J. P. (2012). Multivariate density estimation using dimension reducing information and tail flattening transformations for truncated and censored data. Annals of the Institute of Mathematical Statistics, 48 (1), 167–192. Butt, Z. and Haberman, S. (2004). Application of frailty-based mortality models using Generalized Linear Models. ASTIN Bulletin, 34 (1), 175–197. Cairns, A. J. G., Blake, D. and Dowd, K. (2006). A two factor model for stochastic mortality and parameter uncertainty: theory and calibration. The Journal of Risk and Insurance, 73 (4), 687–718. Cairns, A. J. G., Blake, D., Dowd, K., Coughlan, G. D., Epstein, D., Ong, A. and Balevich, I. (2009). A quantitative comparison of stochastic mortality models using data from England and Wales and the United States. North American Actuarial Journal, 13 (1), 1–35. Cairns, A. J. G., Blake, D., Dowd, K., Coughlan, G. D., Epstein, D. and Khalaf-Allah, M. (2011). Mortality density forecasts: an analysis of six stochastic mortality models. Insurance: Mathematics and Economics, 48, 355–367. Currie, I. D. (2006). Smoothing and forecasting mortality rates with P-splines. Talk given at the Institute of Actuaries, June 2006. http://www.actuaries.org.uk Dowd, K., Cairns, A. J. G., Blake, D., Coughlan, G. D., Epstein, D. and Khalaf-Allah, M. (2010a). Evaluating the goodness of fit of stochastic mortality models. Insurance: Mathematics and Economics, 47, 255–265. Dowd, K., Cairns, A. J. G., Blake, D., Coughlan, G. D., Epstein, D. and Khalaf-Allah, M. (2010b). Backtesting stochastic mortality models: an ex post evaluation of multi-period ahead density forecasts. North American Actuarial Journal, 14 (3), 281–298. Forfar, D. O., McCutcheon, J. J. and Wilkie, A. D. (1988). On graduation by mathematical formula. Journal of the Institute of Actuaries, 115, 1–149. G´ amiz-P´ erez, M. L., Mart´ ınez-Miranda, M. D. and Nielsen, J. P. (2013a). Smoothing survival densities in practice. Computational Statistics and Data Analysis, 58 (1), 368–382. G´ amiz P´ erez, M. L, Mammen, E., Mart´ ınez Miranda, M. D. and Nielsen, J. P. (2013b). Do-validating local linear hazards. Submitted preprint. G´ amiz-P´ erez, M. L., Janys, L., Mart´ ınez-Miranda, M. D. and Nielsen, J. P. (2013c). Smooth marker dependent hazard estimation in praxis. Computational Statistics and Data Analysis,Forthcoming. Guill´ en, M., Nielsen, J. P., Scheike, T. H. and P´ erez-Mar´ ın, A. M. (2012). Time-varying effects in the analysis of customer loyalty: A case study in insurance. Expert Systems with Applications, 39 (3), 3551– 3558. Haberman, S. and Renshaw, A. E. (2011). A comparative study of parametric mortality projection models. Insurance: Mathematics and Economics, 48, 35–55. Horiuchi, S. and Coale, A. J. (1990). Age patterns of mortality for older women: an analysis using the agespecific rate of mortality change with age. Mathematical Population Studies, 2 (4), 245–267. Hougaard, P. (1984). Life table methods for heterogeneous populations: distributions describing the heterogeneity. Biometrika, 71 (1), 75–83. Jarner, S. F. and Kryger, E. M. (2011). Modelling adult mortality in small populations: The SAINT model. ASTIN Bulletin, 41 (2), 377–418. Jones, B. L. (1998). A model for analyzing the impact of selective lapsation on mortality. North American Actuarial Journal, 2 (1), 79–86. Kuang D., Nielsen, B. and Nielsen, J. P. (2011). Forecasting in an extended chain-ladder-type model. Journal of Risk and Insurance, 78 (2), 345–359. Jaap Spreeuw, Jens Perch Nielsen and Søren Fiig Jarner 173 Lee, R. D. and Carter, L. R. (1992). Modeling and forecasting U. S. mortality. Journal of the American Statistical Association, 87 (419), 659–671. Li, J. S.-H., Hardy, M. R. and Tan, K. S. (2009). Uncertainty in mortality forecasting: an extension of the Lee-Carter approach. ASTIN Bulletin, 39 (1), 137–164. Mammen, E., Mart´ ınez-Miranda, M. D., Nielsen, J. P. and Sperlich, S. (2011). Do-validation for kernel density estimation. Journal of the American Statistical Association, 106 (494), 651–660. Mart´ ınez-Miranda, M. D., Nielsen, J. P. and W¨ uthrich, M. V. (2012). Statistical modelling and forecasting of outstanding liabilities in non-life insurance. SORT-Statistics and Operations Research Transactions, 36 (2), 195–218. Nielsen, J. P., Tanggaard, C. and Jones, M. C. (2009). Local linear density estimation for filtered survival data. Statistics, 43 (2), 167–186. Olivieri, A. (2006). Heterogeneity in survival models-applications to pensions and life annuities. Belgian Actuarial Bulletin, 6 (1), 23–39. Pinquet, J., Guill´ en, M. and Ayuso, M. (2011). Commitment and lapse behavior in long-term insurance: a case study. Journal of Risk and Insurance, 78 (4), 983–1002. Plat, R. (2009). On stochastic mortality modelling. Insurance: Mathematics and Economics, 45, 393–404. Renshaw, A. E. and Haberman, S. (2006). A cohort-based extension to the Lee-Carter model for mortality reduction factors. Insurance: Mathematics and Economics, 38, 556–570. Vaupel, J. W., Manton, K. G. and Stallard, E. (1979). The impact of heterogeneity in individual frailty on the dynamics of mortality. Demography, 16, 439–454. Wang, S. and Brown, R. L. (1998). A frailty model for projection of human mortality improvements. Journal of Actuarial Practice, 6, 221–241.