# Parametric and Nonparametric Bootstrap in Actuarial Practice

←

**Page content transcription**

If your browser does not render page correctly, please read the page content below

Parametric and Nonparametric Bootstrap in Actuarial Practice ∗ Krzysztof Ostaszewski and Grzegorz A. Rempala University of Louisville March 20, 2000 Abstract Uncertainty of insurance liabilities has always been the key issue in actuarial theory and practice. This is represented for instance by study and modeling of mortality in life insurance, and loss distributions in traditional actuarial science. These models have evolved from early simple deterministic calculations to more sophisticated probabilistic ones. Such probabilistic models have been traditionally built around parameters characterizing certain probability dis- tributions, e.g., Gompertz’s model of force of mortality, or parametric models of the yield curve process. Modern actuarial science has introduced probabilistic models for all input variables in study- ing insurance firm, and in the whole company models. These new methodologies are based on the theoretical work in mathematical finance which shows that the market, or fair value of insurance liabilities, and indeed, the market value of the insurance firm, can be determined us- ing the general approach developed for contingent claims valuation. While this is theoretically appealing and justified, the central dilemma in modeling insurance company, i.e., its assets and liabilities, is the choice of an appropriate probability distributions, or stochastic processes, governing the evolution of the underlying variables such as interest rates, or asset returns in general, withdrawals or lapses, and mortality. The traditional approaches to this problem have been based on the parametric models. The last two decades have brought about a rich body of new research in nonparametric statistics. This work is intended at showing direct application of a specific nonparametric methodology in actuarial models. The methodology researched here is that of the bootstrap, and more generally, resampling. We develop the bootstrap model alternative on both the asset and liability side. First, we show how bootstrap can be used successfully in enhancing a parametric mortality law suggested by Carriere (1992). Next, we develop a whole company asset-liability model to study a bootstrap alternative to lognormal and stable Paretian models of interest rate process. The results indicate that bootstrap can be instrumental in understanding the rich structure of random variables on the asset and liability sides of an insurance firm balance sheet, and in error estimation for the existing models. 1 Introduction Actuarial science has always been concerned with uncertainty and the process of modeling it. In fact, the uncertainty of insurance liabilities may be considered the very reason for the creation ∗ The authors graciously acknowledge support from the Actuarial Education and Research Fund 1

of actuarial science. This is clearly represented in a study of mortality and other decrements in life, disability, and health insurance, and loss distribution in property/casualty insurance. The models used in such studies began as tabular ones, and then developed into a study of probability distributions. Probability distributions have been traditionally characterized by some numerical parameters. Examples of such probability distributions include Makeham of Gompertz laws of mortality, or the use of gamma distributions in modeling losses, or the omnipresent use of normal approximations (including the lognormal distribution modeling of interest rates). The last two decades have brought about a new vast body of statistical research in nonparametric approaches to modeling uncertainty, in which not the individual parameters of the probability distribution, but rather the entire distribution is sought based on empirical data available. There are varying approaches in nonparametric statistics. For instance, Klugman, Panjer and Willmot (1998) give some examples of the use of kernel density estimation, or nonparametric Bayesian estimation. One of the key new methodologies of nonparametric statistics is the concept of the bootstrap, also knows under a somewhat broader term of resampling. In this work, we present the basic ideas of the bootstrap and resampling, and show how promising its applications in actuarial modeling can be. The paper is organized as follows. In the next section we give a brief overview of the basic ideas of the bootstrap methodology in the context of parametric and non-parametric estimation as well as time series-based inferences. The subsequent two sections of the paper are devoted to the somewhat detailed discussion of two particular examples, arising from the actuarial modeling problems. In Section 3 we we look at mortality estimation, studying a mortality law introduced by Carriere (1992), generalizing on the classical Gompertz mortality law. In Section 4, we proceed to a whole company asset-liability model, comparing cash flow testing results for an annuity company under two parametric assumptions of interest rate process and the method based on the nonparametric, bootstrap-based approach to analyzing time series data. Section 5 contains some discussion of the obtained results whereas Section ?? offers some concluding remarks. 2 The Bootstrap Before we discuss the applications of resampling methodology, let us first give a brief overview of its basic ideas. 2.1 The Concept The concept of the bootstrap was first introduced in the seminal piece of Efron (1979) as an attempt to give some new perspective to an old and established statistical procedure known as jackknifing. Unlike jackknifing which is mostly concerned with calculating standard errors of statistics of interest, Efron’s bootstrap has been set to achieve more ambitious goal of estimating not only the standard error but also the distribution of a statistic. In his paper Efron has considered two types of bootstrap procedures useful, respectively, for nonparametric and parametric inference. The nonparametric bootstrap relies on the consideration of the discrete empirical distribution generated by a random sample of size n from an unknown distribution F . This empirical distribution Fbn assigns equal probability to each sample item. In the parametric bootstrap setting, we consider F to be a member of some prescribed parametric family and obtain Fbn by estimating the family parameter(s) from the data. In each case, by generating an iid random sequence, called a resample or pseudo-sequence, from the distribution Fbn , or its appropriately smoothed version, we can arrive at new estimates of various 2

parameters or nonparametric characteristics of the original distribution F . This simple idea is at the very root of the bootstrap methodology. In 1981 Bickel and Freedman (1981) formulated conditions for consistency of the bootstrap. This resulted in further extensions of the Efron’s methodology to a broad range of standard applications, including quantile processes, multiple regression and stratified sampling. They also argued that the use of bootstrap did not require theoretical derivations such as function derivatives, influence functions, asymptotic variances, the Edgeworth expansion, etc. Singh (1981) made a further point that the bootstrap estimator of the sampling distribution of a given statistic may be more accurate than the traditional normal approximation. In fact, it turns out that for many commonly used statistics the bootstrap is asymptotically equivalent to the one- term Edgeworth expansion estimator, usually having the same convergence rate, which is faster than the normal approximation. In many recent statistical texts, the bootstrap is recommended for estimating sampling distribu- tions, finding standard errors, and confidence sets. Although the bootstrap methods can be applied to both parametric and non-parametric models, most of the published research in the area has been concerned with the non-parametric case since that is where the most immediate practical gains might be expected. Bootstrap procedures are available in most recent editions of some popular statistical software packages, including SAS and S-Plus. Recently, several excellent books on the subject of bootstrap, resampling and related techniques have become available. Readers interested in gaining some basic background in resampling are referred to Efron and Tibisharani (1993). For a more mathematically advanced treatment of the subject, we recommend Shao and Tu (1995) which also contains a chapter on bootstrapping time series and other dependent data sets. 2.2 Bootstrap Standard Error and Bias Estimates Arguably, one of the most important practical applications of the bootstrap is in providing concep- tually simple estimates of the standard error and bias for a statistic of interest. Let θ̂n be a statistic based on the observed sample, arriving from some unknown distribution function F . Assume that θ̂n is to estimate some (real-valued) parameter of interest θ, and let us denote its standard error and bias by seF (θ̂n ) and biasF (θ̂n ). Since the form of the statistic θ̂n may be vary complicated the exact formulas for the corresponding bootstrap estimates of standard error (BESE) and bias (BEB) may be quite difficult,if not impossible, to derive. Therefore, one usually approximates both these quantities with the help of the multiple resamples. The approximation to the bootstrap estimate of standard error of θ̂n suggested by Efron (1979) is given by ( B )1/2 X ∗ bB = se [θ̂ (b) − θ̂n∗ (·)]2 /(B − 1) (1) b=1 where θ̂n∗ (b) is the original statistic θ̂n calculated from the b-th resample (b = 1, . . . , B), θ̂n∗ (·) = PB ∗ b=1 θ̂n (b)/B, and B is the total number of resamples (each of size n) collected with replacement from the empirical estimate of F (in parametric or non-parametric setting), By the law of large numbers b B = BESE(θ̂n ), lim se B→∞ and for sufficiently large n we expect BESE(θ̂n ) ≈ seF (θ̂n ). d B based on B resamples Similarly, for BEB one can use its approximation bias 3

X B dB = bias θ̂n∗ (b)/B − θ̂n (2) b=1 Let us note that B, total number of resamples, may be taken as large as we wish, since we are in complete control of the resampling process. For instance, it has been shown that for estimating BESE, B equal to about 250 typically gives already a satisfactory approximation, whereas for BEB this number may have to be significantly increased in order to reach the desired accuracy (see Efron and Tibisharani 1993 for a discussion of these issues). 2.3 Bootstrap Confidence Intervals Let us now turn into the problem of using the bootstrap methodology to construct confidence intervals. This area has been a major focus of theoretical work on the bootstrap and in fact the procedure described below is by far not the most efficient one and can be significantly improved in both rate of convergence and accuracy. It is, however, intuitively obvious and easy to justify and seems to be working well enough for the cases considered here. For complete review of available approaches to bootstrap confidence intervals for iid data, see Efron and Tibisharani (1992). As in the case of standard error and bias estimates the methodology for bootstrap interval estimation applies as well to time dependent observations as long as we modify our resampling procedure to ether block or circular bootstrap sampling (see Section 2.4 below). Let us consider θ̂n∗ , a bootstrap estimate of θ based on a resample of size n (or, in the case of dependent observations, on k blocks of length l, see Section 2.4 below) from the original data sequence X1 , . . . , Xn , and let G∗ be its distribution function given the observed series values G∗ (x) = P rob{θ̂n∗ ≤ x|X1 = x1 , . . . , Xn = xn }. (3) Recall that for any distribution function F and τ ∈ (0, 1) we define the τ -th quantile of F (sometimes also called τ -th percentile) as F −1 (τ ) = inf{x : F (x) ≥ τ }. The bootstrap percentiles method gives G−1 −1 ∗ (α) and G∗ (1 − α) as, respectively, lower and upper bounds for 1 − 2α confidence interval for θ̂n . Let us note that for most statistics θ̂n the form of a distribution function of the bootstrap estimator θ̂n∗ is not available. In practice, G−1 −1 ∗ (α) and G∗ (1 − α) are approximated by generating B pseudo-sequences (X1 , . . . , Xn ), calculating the corresponding values of θ̂n∗ (b) for b = 1, . . . , B, ∗ ∗ and then finding the empirical percentiles. In this case the number of resamples B usually needs to be quite large; in most cases it is recommended that B ≥ 1000. 2.4 Dependent Data It is not difficult to show that the Efron’s bootstrap procedure fails when the observed sample points are not independent. The extension of the bootstrap method to the case of dependent data was first considered by Künsch (1989) who suggested a moving block bootstrap procedure which takes into account the dependence structure of the data by resampling blocks of adjacent observations rather than individual data points. The method was shown to work reasonably well but suffered a drawback related to the fact that for fixed block and sample sizes observations in the middle part of the series had typically a greater chance of being selected into a resample than the observations close to the ends. This happens because the first or last block would be often shorter. To remedy this deficiency Politis and Romano (1992) suggested a method based on circular blocks, i.e. on 4

wrapping the observed time series values around the circle and then generating consecutive blocks of bootstrap data from the circle. In the case of the sample mean this method known as a circular bootstrap again was shown to accomplish the Edgeworth correction for dependent data. Let X1 , . . . , Xn be a (strictly) stationary time series, and as before let θ̂n = θ̂n (X1 , . . . , Xn ) denote a real-valued statistic estimating some unknown parameter interest θ. Given the observations X1 , . . . , Xn and an integer l (1 ≤ l ≤ n) we form the l-blocks Bt = (Xt , . . . , Xt+l−1 ) for t = 1, . . . , n − l + 1. In moving blocks procedure of Künsch (1989) the bootstrap data (pseudo-time series) (X1∗ , . . . , Xn∗ ) is obtained by generating the resample of k blocks (k = [n/l]) and then considering all individual (inside-block) pseudo-values of the resulting sequence B1∗ . . . , Bk∗ . The circular bootstrap approach is somewhat similar except that here we generate a subsample ∗ ∗ B̃1 . . . , B˜k from the l-blocks B̃t = (X̃t , . . . , X̃t+l−1 ) for t = 1, . . . , n where ( Xt if t = 1, . . . , n, X̃t = Xt−n t = n + 1, . . . , n + l − 1. In either case, the bootstrap version of θ̂n is θ̂n∗ =θ̂n (X1∗ , . . . , Xn∗ ). The idea behind both the block and circular bootstrap is simple: by resampling blocks rather than original observations we preserve the original, short-term dependence structure between the observations although not necessarily the long-term one. For the data series for which the long term dependence is asymptotically negligible is some sense (e.g., α-mixing sequences with appropriately fast mixing rate) the method will produce consistent estimators of the sampling distribution of θ̂n −θ as long as l, n → ∞ at the appropriate rate (see e.g., Shao and Tu Chapter 9 for the discussion of the appropriate assumptions). In particular, it follows that the formulas for the approximate BESE and BEB given by (1) and (2), respectively as well as the method of percentiles for confidence estimation are still valid with the block and circular bootstrap. Let us also note that for l = 1 both methods reduce to the Efron’s bootstrap procedure for iid random variables. In the next part of the paper we demonstrate the use of the above bootstrap techniques in both parametric and non-parametric setting. 3 Modeling US Mortality Tables In this section we consider an application of the bootstrap techniques to estimating the variability of the parameters of a general mortality law proposed by Carriere (1992). 3.1 Carriere Mortality Law The research into a law of mortality has had a long history going back the famous formula of De Moivre. Probably the most popular parametric model of mortality (proposed by Benjamin Gompertz) has been the one that models the force of mortality µx = Bcx (4) as an exponential function. Whereas Gompertz’s model typically fits observed mortality rates fairly well at the adult ages it exhibits some deficiency in modeling the early ones. This is mostly due to the fact that in many populations the exponential growth in mortality at adult ages is preceded by a sharp mortality fall in early childhood and the hump at about age 23 (cf. e.g., Tenenbein and Vanderhoof 1980). This early years behavior of the mortality force is not easily modeled with the 5

formula (4). There has been many attempts in the actuarial literature to modify Gompertz’s law in order to obtain a better approximation for the early mortality patterns. It has been observed in Heligman and Pollard (1980) that by adding additional components to Gompertz model one can obtain a mortality law that fits early mortality patterns much better. This idea has been further extended by Carriere 1992 [referred to in the sequel as CR] who has suggested using a mixture of the extreme-value distributions: Gompertz, Inverse Gompertz, Weibull and Inverse Weibull, to model the mortality curve of a population. Namely, he suggested a mortality law based on a survival function s(x) given as the mixture X 4 s(x) = ψi si (x) (5) i=1 P4 where 0 ≤ ψi for i = 1 . . . , 4, i=1 ψi = 1 are the mixture coefficients and si (x) for i = 1, . . . 4 are survival functions based on the Gompertz, Inverse Gompertz, Weibull, and Inverse Weibull distribution functions, respectively, with the corresponding location and scale parameters mi and σi . More precisely, s1 (x) = exp(e−m1 /σ1 − e(x−m1 )/σ1 ) (6a) −(x−m2 )/σ2 1 − exp(−e ) s2 (x) = (6b) 1 − exp(−em 2 /σ 2 ) s3 (x) = exp(−(x/m3 )m3 /σ3 ) (6c) −m4 /σ4 s4 (x) =1 − exp(−(x/m4 ) ) (6d) In Carriere mortality model the Gompertz (6a) and Inverse Gompertz (6b) curves, which both are derived from (4), are used to model the mortality of adulthood, whereas the Weibull and Inverse Weibull curves are used to model the mortality of childhood and early adolescence. The introduction of these last two survival functions allows in particular to fit correctly the early years decrease in mortality as well as the hump around age 20, which are both present in many populations mortality patterns. As it has been demonstrated empirically in CR the model (5) fits quite well the patterns of mortality for the entire US population as well as separately US male and female mortality laws. Due to the interpretation of the survival components of the curve (5) Carriere model, unlike that of Heligman and Pollard (1980), gives also the estimates of childhood, adolescence and adulthood mortality by means of the mixture coefficients ψi . In the next section we discuss a way of arriving at these estimates. 3.2 Fitting the Mortality Curve As an example, we fit the general mortality model (5) to the Aggregate 1990 US Male Lives Mortality Table (US Health Dept. 1990). This particular table is of special actuarial interest because it was a data source in developing 1994 GAM and 1994 UP tables. In order to fit the model (5) to the observed mortality pattern, some fitting criterion (a loss function) is needed. In his paper CR has considered several different loss functions but here for the sake of simplicity we will restrict ourselves to only one of them given by 100 µ X ¶2 q̂x 1− , (7) x=0 qx 6

Survival component Location par (mi ) Scale par (σi ) Mixture par (ψi ) s1 (·) 80.73 11.26 0.944 s2 (·) 42.60 14.44 0.020 s3 (·) 0.30 1.30 0.013 s4 (·) 23.65 7.91 0.022 Table 1: Estimates P3 of the Carriere survival model parameters for 1990 US Male Mortality Tables. ψ4 = 1 − i=1 ψi where qx is a probability of death within a year for a life aged x, as obtained from the US table and q̂x is a corresponding probability calculated using the survival function (5). The ages above x = 100 were not included in (7) since, similarly to CR, we have found that the mortality pattern for the late ages in the table had been based on Medicare data and hence might be not representative to the mortality experience of the entire US male population. The parameters of the survival functions (6) as well as the corresponding mixing coefficients ψi ’s were calculated by minimizing the loss function (7). The calculations were performed with a help of the SOLVER add-on function in the Microsoft Office 97 Excel software. The estimated values of the parameters mi , σi2 , and ψi for i = 1, . . . 4 are given in the Table 1. As we can see from the table the modes mi of the survival components in the Gompertz and the Inverse Gompertz case are equal to 80.73 and 42.60 respectively, whereas in the Weibull and the Inverse Weibull case they are equal to 0.3 and 23.65, respectively. This information reveals that, as intended, both Gompertz components model the later age mortality, and Weibull as well as Inverse Weibull components model the mortality of childhood and adolescence. Moreover, since ψ1 = 0.94 it follows that most of the deaths in the considered population are due to the Gompertz component. A plot of the fitted curve of ln(q̂x ) using the estimated values of the parameters given in Table 1 along with the values of ln(qx ) for x = 0, . . . , 100 is presented in Figure 1. As we can see from the plot the fit seems to quite good except perhaps between the ages of 6 and 12. Let us also note the remarkable resemblance of our plot in Figure 1 with that presented in Figure 4 of CR, which was based on the 1980 US Population Mortality Table. 3.3 Statistical properties of the parameter estimates in Carriere mortal- ity model The statistical properties of estimators of the unknown parameters ψi i = 1 . . . 3 and mi , σi i = 1 . . . , 4 obtained by minimizing the function (7) (we refer to them in the sequel as the Carriere or CR estimators) are not immediately obvious. For instance, it is not clear whether or not they are consistent and what is their asymptotic efficiency. Since these questions are related to the problem of consistency of the bootstrap estimation, we will discuss them briefly here. In order to put our discussion in the context of the statistical inference theory we need to note an important difference between the loss function (7) and the usual, weighted least-square loss. Namely, in our setting the US mortality rates qx for x = 1, . . . , 100 are not independent observations but rather are themselves estimated from the large number n of the deaths observed in the population.1 Since the qx ’s are estimated based on the same number n of the observed deaths, they are negatively correlated and hence the loss function (7) treated as a random function, 1 Examining the methodology of creating the US table used as an example in this section we found that in fact the qx ’s were obtained by interpolation with the help of Beer’s 4-th degree oscillatory formula. 7

Figure 1: A plot of ln(qx ) from US Life Table (dashed line) and of ln(q̂x ) using formulas (6) with the parameters values given in Table 1. is not a sum of independent identically distributed random variables. In view of the above, the usual statistical methodology of M -estimation does not apply and in order to obtain asymptotics for the CR estimators we need to develop a different approach which is briefly sketched below. Under the assumption that the true mortality curve of the population indeed follows model (5) with some unknown parameters ψi0 , (i = 1, . . . , 3) and m0i , σi0 (i = 1, . . . 4) one may consider CR estimators associated with the i-th survival component (i = 1, . . . , 4) as functions of the vector of mortality rates qx for x = 1, . . . , 100 and expand them around the true parameters ψi0 , m0i , σi0 using the multivariate Taylor’s Theorem. The expansion is valid in view of the Implicit Function Theorem which guarantees that all the local minimizers of (7) are locally continuously differentiable functions of the qx ’s as long as the first derivatives of q̂x = q̂x (mi , σi , ψi ) are locally continuous and do not vanish (that is, the appropriate Jacobian function is no equal to zero, at least for one age x) around the point ψi0 , m0i , σi0 . This method reveals that CR estimates are consistent and that their join distribution is asymptotically multivariate normal with some covariance matrix of the entries depending upon the values of the true parameters ψi0 , (i = 1, . . . , 3), m0i , σi0 (i = 1, . . . 4) and two first derivatives of the q̂x ’s with respect to ψi , mi , σi . 3.4 Assessment of the Model Accuracy with Parametric Bootstrap One way of assessing the variability of CR estimates would be to use the normal approximation outlined in the previous section. However, since the asymptotic variance of this approximation depends upon the unknown parameters which would have to be estimated from the data, the derivation of the exact formula for the covariance matrix as well as the assessment of the performance of variance estimators would be required. In addition, the normal approximation would not account 8

for the possible skewness of the marginal distributions of ψi , mi , and σi .. Therefore, as an alternative to this traditional approach, we propose to apply the bootstrap method instead. Since the bootstrap typically corrects for the skewness effect we would hope that it might provide in our setting a better approximations than the methods based on the normal theory. Let us note that the methods of the nonparametric bootstrap cannot be applied here since we do not know the empirical distribution Fbn of the observed deaths. However, we do know the estimated values of the parameters in the model (5), since these are obtained by minimizing the expression (7) which is concerned only with the quantities qx . Hence, assuming that the true population death distribution follows the model (5) we may use the values given in Table 1 to obtain our resampling distribution Fbn . Once we have identified the resampling distribution, the bootstrap algorithm for obtaining approximate variances, biases, and confidence intervals for CR estimators is quite simple. We use a computer to generate a large (in our example 20,000 points) random sample from Fbn and apply interpolation2 to obtain the sequence of the pseudo mortality rates qx∗ for x = 1, . . . , 100. In order to find the corresponding bootstrap replication ψi∗ i = 1 . . . 3 and m∗i , σi∗ i = 1 . . . , 4 we minimize the loss function (7), where now the qx ’s are replaced by the qx∗ ’s. The above procedure is repeated a large number (B) times and the resulting sequence of bootstrap replications ψi∗ (b), m∗i (b), σi∗ (b) for b = 1, . . . , B is recorded. The plot of a single replication of the log of qx∗ ’s curve and the corresponding log of q̂x∗ ’s curve, i.e., log of the fitted survival curve (5) with the parameters ψi∗ , m∗i , σi∗ is presented in Figure 2. In order to find the corresponding bootstrap estimates of standard error and bias for each one of the 11 CR estimators in Table 1 we use the formulas (1) and (2). The corresponding bootstrap confidence interval are calculated as well, with the help of the bootstrap percentile method. The general appropriateness of the above outlined bootstrap procedure is not immediately ob- vious and in fact requires some theoretical justification as it doe not follow immediately from the general theory. Since the mathematical arguments needed to do so get quite technical and are beyond the scope of the present work, let us just only briefly comment that, generally speaking, for the mere consistency of our bootstrap algorithm, the argument similar to that used to argue the consistency of CR estimators may be mimicked. However, in order to show that our bootstrap estimators in fact have an edge over the normal approximation method a more refined argument based on the Edgeworth expansion theory is needed. For some insight into this and related issues see, for instance, Hall (1992). The final results based on our bootstrap approximations to the standard error, bias and confi- dence bounds for all 11 CR estimators given in Table 1 are presented in Table 2. As we can see from these results it seems that only the estimates of the Gompertz component (m1 , σ1 , ψ1 ) of the mortality model enjoy a reasonable degree of precision in the sense that their standard errors and bias are small (relative to the corresponding numerical value), and their confidence intervals short. For all the other component both the corresponding values for standard errors and the lengths of the confidence interval, indicate their apparent lack of precision. In particular, the lower bound of the mixture coefficient ψ2 corresponding to the Inverse Gompertz component (6b) equals zero, which could indicate that the Inverse Gompertz component is simply not present in the considered mortality pattern. In fact, CR fitting the model (5) to US 1980 Mortality data, which has a similar pattern to the one we have considered here, did take ψ2 = 0. Hence, by providing the 95% CI’s for 2 In order to maximize the efficiency of the procedure one should use here the same method of interpolation as the one used to create original values of the qx ’s in the table, however, in our example we simply applied linear interpolation to obtain the values of the pseudo-empirical distribution at the integer ages and then to calculate the qx∗ ’s 9

Figure 2: A plot of ln(qx∗ ) (rough solid line) from a single 20,000 point resample of a parametric model (dashed line) and of ln(q̂x∗ ) latest marked by a smooth solid curve). The dashed curve is the original resampling curve with its parameters values given in Table 1. the estimates of the ψi ’s our bootstrap method also suggests the general form of the mixture in the Carriere law.3 For illustration purposes the distribution of B = 1000 bootstrap replications for the parameter m1 along with its 95% confidence interval is presented in Figure 3. Let us note that the histogram of the replication values seems to follow the normal curve (although it is somewhat skewed to the left), which is consistent with our theoretical derivations. As we have seen above the parametric bootstrap algorithm can be helpful in obtaining a more complete information about the fitted mortality model and its parameters. The drawback of the algorithm presented above is mainly in its very high computational cost, since for every resample of the sequence of 11 values of the parameters we need to first obtain a resample of the fitted mortality rates q̂x , x = 1 . . . , 100. It would be interesting to develop some better methods of bootstrapping in the problems involving the loss functions of type (7), either with the help of the residuals bootstrap method (see Efron and Tibisharani 1993) or maybe some techniques aiming at increasing the efficient of resampling, like e.g. importance method (see Shao and Tu 1995). 3 Although in this paper we are concerned mostly with the point estimation issues and do not discuss the hypothesis testing, it is not difficult to see that the nonparametric bootstrap method considered here could be also used to perform a formal test of significance of ψi = 0 (in a similar fashion to the usual test of the coefficients in the regression models). For more on the topic of hypothesis testing with the bootstrap, see e.g., Efron and Tibisharani (1993). 10

Figure 3: the distribution of B = 1000 bootstrap replications for the parameter m1 along with its 95% confidence interval. The vertical line in the center shown the actual fitted value of m1 as given in Table 1. 4 Cash flow testing company model As our second example of the use of the resampling methods, we discuss the applications of the bootstrap to modeling a company cash flow as a function of the interest rates process. Modern financial methodology of pricing contingent claims, began with the seminal work of Black and Scholes (1973), has had increasing influence on the actuarial practice. The emerging consensus in the actuarial profession calls for consistent use of cash flow testing in analyzing solvency, profitability, fair value of liabilities, and the value of an insurance firm (see, e.g., Vanderhoof and Altman, 1998, also Ostaszewski, 2000). One of the key features of such methodologies is the presence of numerous random variables and stochastic processes in the models, in addition to the standard assumption of randomness of losses or mortality. The most important random effect missing in traditional methodology has been interest rates, or capital asset returns in general. One could also point out randomness of other factors, mainly responses of consumers to economic incentives (lapse rates, growth of new business, mortgage prepayment rates, etc.). However, the effect of interest rates is so pronounced, that the 1991 Standard Valuation Law, and common actuarial practice, now include some form of cash flow testing under varying interest rates scenarios as a standard. 4.1 Interest Rates Process While there exists a consensus as to the necessity of using interest rates scenarios, and to account for some form of their randomness, one can hardly argue that such a consensus is formed concerning the actual stochastic process governing the interest rates process. Becker (1991) analyzes one common 11

m1 σ1 ψ1 m2 σ2 ψ2 Value (from Table 1) 80.73 11.26 0.944 42.60 14.44 0.020 SE 0.245 0.278 0.018 8.257 4.788 0.015 Bias 0.096 -0.120 -0.003 -1.266 -1.653 0.005 Lower Bound 95% CI 80.414 10.521 0.900 26.222 3.007 0.000 Upper Bound 95% CI 81.336 11.608 0.965 54.531 19.932 0.056 m3 σ3 ψ3 m4 σ4 ψ4 Value (from Table 1) 0.30 1.30 0.013 23.65 7.91 0.022 SE 0.303 2.268 0.001 3.011 3.122 N/A Bias 0.037 0.254 0.000 -0.179 -0.253 N/A Lower Bound 95% CI 0.053 0.359 0.011 17.769 3.048 N/A Upper Bound 95% CI 1.290 10.019 0.015 30.000 13.504 N/A Table 2: Bootstrap estimates of the standard error, bias and 95% CI for the estimators of the Carriere survival model parameters in 1990P US Male Mortality Tables presented in Table 1. The 3 mixture coefficient ψ4 is given by ψ4 = 1 − i=1 ψi and hence is not a model parameter. model derived from the process dIt = It µdt + It σdW (8) with It denoting the interest rate prevailing at time t, W denoting the standard Wiener process, with µ and σ being the drift and volatility parameters. Becker tests the hypothesis, implied by substituting µ = 0 and σ = const, that the distribution of µ ¶ It+1 Jt = ln It is normal with mean zero and constant variance. This means that the ratio of interest rates in two consecutive periods has the lognormal probability distribution. However, Becker (1991) rejects the hypothesis of normality and even independence of the Jt ’s based on the empirical evidence. Of course, (1) is not the only stochastic differential equation proposed to describe interest rate process. Many other processes have been put forth in this very active area of modern mathematical finance. Hull (1997) provides an overview of various stochastic differential equations studied in this context. De La Grandville provides a straightforward argument for approximate lognormality of the variable (1 + rate of return), based on the Central Limit Theorem. But Klein (1993) discusses various reasonings underlying lognormality assumption, and points out that mutual independence of con- secutive returns, finite variance of one period return distribution, invariance under addition, and the functional form of the distribution, are very convenient but often unreasonable assumptions. Indeed, many of them have been rejected by a series of recent nonparametric studies of capital markets, discussed by Cerrito, Olson and Ostaszewski (1998). While Klein (1993) argues convinc- ingly against finite variance assumption, the work of Cerrito, Olson and Ostaszewski (1998), among others (notably Fama and French, 1988), aims at the evidence against independence of returns in consecutive periods. Both take away key assumptions of the Central Limit Theorem. Klein (1993) [referenced as KL in the sequel] proposes to replace lognormal distribution by the Stable Paretian, and provides a comprehensive comparison of whole company models under these two hypotheses about interest rates process. As these differing assumptions produce drastically different results 12

based on the same data, KL argues that a valuation actuary must give proper consideration to the possibility that the stochastic dynamics of the interest rate process differ significantly from the simple assumption of lognormality. 4.2 Modeling Interest Rates with Nonparametric Bootstrap of Depen- dent Data. Let us note that, unlike in the case of the US Mortality Tables, the historical, empirical realizations of the interest rates processes are typically available, which allows us to use the non-parametric bootstrap methods. In the following sections of this paper we investigate a fully non-parametric model of a company cash flow as a function of the interest rates process {Jt }. In particular, we do not assume independence of the Jt ’s, relaxing that assumption to only have interest rates process stationarity. Furthermore, we propose that it may be unreasonable to assume one functional form of the probability distribution. The choice of the non-parametric model here is not merely for the sake of example. In our opinion, this allows to proceed to deeper issues concerning the actual forms of uncertainty underlying the interest rate process, and other variables studied in dynamic financial analysis of an insurance firm. In our view, one cannot justify fitting convenient distributions to data and expect to easily survive the next significant change in the marketplace. What works in practice, but not in theory, may be merely an illusion of applicability provided by powerful tools of modern technology. If one cannot provide a justification for the use of a specific parametric distribution, than a nonparametric alternative should be studied, at least for the purpose of understanding firm’s exposures. In this work, we will study the parametric and nonparametric bootstrap methodologies as applied to an integrated company model of KL. 4.3 Assumptions for the Cash-Flow Analysis While there are many variables involved in a company model, in our example we will limit ourselves to considering a nonparametric estimate of the distribution of 10-th year surplus value as a stochastic function of the underlying interest rates process. We will also assume a simple one-parameter structure of the yield curve. Further work will be needed for extension of this idea to nonparametric estimates of the entire yield curve, and the liability side variables. Stanton (1997) provides one possible nonparametric model of term structure dynamics and the market price of interest rate risk. In addition to the simple one parameter model of the yield curve, we also make a simplifying assumption of parallel shifts of all interest rates derived from the basic rate (long term Treasury rate). This follows the design of the model company in KL. Clearly, actuarial practice will eventually require greater sophistication of this model. From our perspective, this requires extending bootstrap techniques to several dependent time series (e.g., interest rate series at various maturities). We believe this is definitely an area for important further research. The model in this work, following that of KL, also assumes a parametric functional form for lapses, withdrawals, and mortality, as well as for prepayments for mortgage securities in the asset portfolio. Again, we believe that further research into nonparametric alternatives to lapse, mor- tality, and prepayment models, would be very beneficial to our understanding of the probability distribution structure of those phenomena. It should be noted that any research in these areas will require extensive data, which, unlike the interest rates and capital asset returns, is not always easily available. Overall, we would like to stress that our work in this area has to be, due to numerous limitations, limited, but we do believe that the actuarial profession would greatly benefit from formulation of 13

nonparametric, bootstrap-based models of whole company, incorporating models for all types of uncertainty, i.e., asset returns, lapses, mortality, prepayments, and other random phenomena. 4.4 Company Model Assumptions Following the design of KL we create a model company for cash flow analysis comparison of the traditional lognormal model of interest rates, and stable Paretian distribution discussed by KL, and our nonparametric model based on resampling. This will allow us to have a direct comparison of parametric and nonparametric models, and make conclusions based on such comparisons. The company is studied effective December 30, 1990. It offers a single product: deferred annuity. The liability characteristics are as follows: • Number of policies: 1,000 • Fund Value of Each Policy: 10,000 • Total Reserve: 10,000,000 • Surrender Charge: None • Minimum Interest Rate Guarantee: 4% • Reserve Method: Reserve Equals to Account Balance The assets backing this product are: • 8,000,000 30-year 9.5% GNMA mortgage pools • 2,000,000 1 Year Treasury Bills for a total asset balance of 10,000,000 effectively making initial surplus equal to zero. Again matching KL analysis, we will assume the following initial yields in effect as of December 31, 1990: • 1 Year Treasury Bills: 7.00% • 5 Year Treasury Notes: 7.50% • 30 Year Treasury Bonds: 8.25% • Current coupon GNMAs: 9.50% Treasury yields are stated on a bond-equivalent basis, while GNMAs yields are nominal rates, compounded monthly. The following interest rates are given on 5 Year Treasury Notes at indicated dates: • December 31, 1989: 7.75% • December 31, 1988: 9.09% • December 31, 1987: 8.45% • December 31, 1986: 6.67% 14

KL determines the parameters of the lognormal and Stable Paretian distribution for 30 year Trea- suries yields based on the historical yields for the period 1953-1976 and 1977-1990. He proceeds then to generate 100 sample 30 year Treasuries yields interest rates scenarios for 10 year period. As in KL, the following assumptions about the model company were made: • Lapse rate formula (based on competition and credited interest rates, expressed as percent- ages) is given by ( 2 (w) 0.05 + 0.05 [100(icomp − icred )] if icomp − icred > 0; qt = 0.05 otherwise; with an overall maximum of 0.5. • Credited interest rate (icred ) is a currently anticipated portfolio yield rate for the coming year on a book basis less 150 basis points. In addition to always crediting no less than the minimum interest rate guarantee, the company will always credit a minimum of icomp − 2%. • Competition interest rate (icomp )is the greater of 1 Year Treasury Bill nominal yield to ma- turity and the rate 50 basis points below 5 year rolling average of 5 Year Treasury Bond’s nominal yield to maturity. • Expenses and Taxes are ignored. • There is no annuitization. • Mortgage prepayment rate (ratet ) for the 30 year GNMA pools, calculated separately for each year’s purchases of GNMAs: ( 2 0.05 + 0.03 [100(icoup,t − icurr )] + 0.02 [100(icoup,t − icurr )] if icomp − icred > 0; ratet = 0.05 otherwise; with an overall maximum of 0.40. Here, icoup,t denotes the coupon on newly issued 30 year GNMAs, issued in year t, and icurr denotes the coupon rate on currently issued 30 year GNMAs. which is assumed to shift in parallel with the yields on 30 year Treasury Bonds. We assume the same, i.e., parallel shifts, about the 1 Year Treasury Bill rates, and the 5 Year Treasury Bond rates. • Reinvestment strategy is as follows. If the net cash flow in a given year is positive, any loans are paid off first, then 1-year Treasury Bills are purchased until their book value is equal to 20book value of assets, and then newly issued current coupon GNMAs are purchased. • Loans, if needed, can be obtained at the rate equal to current 1 Year Treasury Bills yield. • Book Value of GNMAs is equal to the present value of the future principal and interest payments, discounted at the coupon interest rate at which they were purchased, assuming no prepayments. • Market Value of GNMAs is equal to the present value of the future principal and interest payments, discounted at the current coupon interest rate, and using the prepayment schedule based on the same rate, assuming that the current coupon rate will be in effect for all future periods as well. 15

• Market Value of Annuities is equal to the account value. All cash flows are assumed to occur at the end of the year. As we had indicated, the projection period is 10 years. We run 100 scenarios and compare them to scenarios run by KL. At the end of 10 years, market values of assets and liabilities are calculated. The distribution of final surplus is then compared. 4.5 Interest Rates Process Assumptions In order to make a comparison of our non-parametric cash-flow analysis with that presented in KL, we have considered the same set of yield rates on 30-years Treasury bonds from years 1977–1990, namely the average annualized yield to maturity on long-term treasury bonds from years 1977-1990. The data is presented in Table 3 below. In the first column of Table 3 we give the average yield to maturity on 30-Years Treasury Bonds. These semiannual (bond-equivalent) rates are then converted into effective annual rates It , for t = 1, . . . 168 where t is the number of months beyond December 1976. The realizations of these converted rates are presented in column 2. Finally, the realized values of the random time series Jt , calculated according to the formula Jt = ln(It+1 /It ) are given in the third column. Let us first note that traditional actuarial assumption that Jt are iid normal do not seem to be satisfied in this case. As shown in Figure 4, the plot of Jt empirical quantiles versus that of a normal random variable indicates that the Jt ’s possibly comes a distribution other than normal. In his approach KL addressed this problem by considering a different, more general form of a marginal distribution for Jt , namely that of a Stable Paretian family, but he still assumed independence of the Jt ’s. However, as it was pointed out in the discussion of his paper by Beda Chan there is clear statistical evidence that the Jt ’s are not, in fact, independent. The simplest illustration of that can be obtained by applying the runs test (as offered for instance by Minitab 8.0 statistical package). There are also more sophisticated general tests for a random walk in a time series data under ARMA(p, q) model as offered by SAS procedure ARIMA (SAS Institute 1996). The results of these tests conducted by us for the data in Table 3 give strong evidence against the hypothesis of independence (both P -values < .001). In view of the above, it seems to be of interest to consider the validity of the stationarity assumption for Jt . Under ARMA(p, q) model again the SAS procedure ARIMA may be used to perform a formal test of non-stationarity either by means of Phillips-Perron or Dickey-Fuller tests. However, a simpler graphical method based on a plot of an autocorrelation function (ACF) and values of the appropriate t-statistics (cf. e.g., Bowerman and O’Connell 1987) may also be used (if we assume that the marginal distribution of Jt has a finite variance). In our analysis we have used both formal testing and graphical methods and found no statistical evidence of non-stationarity of Jt (in conducted tests all P-values < .001). The plot of ACF for Jt is given in Figure 2. In view of the above results, for the purpose of our non-parametric cash-flow analysis we have made the following assumptions about the time series Jt . (A1) It is a (strictly) stationary time series and (A2) has the m-dependence structure, i.e., t-th and s + t + m-th elements in the series are inde- pendent for s, t = 1, 2, . . . and some fixed integer m ≥ 1. The assumption (A1) implies in particular that all the Jt ’s have the same marginal distribution. Let us note, however, that we have not assumed here anything about the form of this marginal distribution and, in particular, we have made no assumptions regarding the existence of any of its moments. Thus Jt could possibly have an infinite variance or even infinite expectation (which is the case for some distributions from the Stable Paretian family). As far as the assumption (A2) 16

17

Table 3: The average yield to maturity on 30-Years Treasury Bonds along with the annualized yield (It ) and the natural logarithm of the ratios of consecutive annualized yield rates (Jt ) 0.13 0.08 0.03 Jt -0.02 -0.07 -0.12 -3 -2 -1 0 1 2 3 Normal Distribution Figure 4: The quantile-to-quantile plot of Jt values given in Table 3 vs the normal distribution is concerned, the ACF plot in Figure 2 indicates that in the case of long-term Treasury Bonds from years 1953-1976, an appropriate value of m is at least 50. The assumption of m-dependence for Jt seems to be somewhat arbitrary, and in fact from a theoretical viewpoint in order for our bootstrap methodology to apply, the condition (A2) may be indeed substantially weakened by imposing some much more general but more complicated technical conditions on the dependence structure of Jt . These conditions in essence require only that Jt and Jt+s are “almost” independent as s increase to infinity at appropriate rate and are known in the literature as mixing conditions. Several types of mixing conditions under which our method of analysis is valid are given by Shao and Tue (1995 p.410) For the sake of simplicity and clarity of presentation, however, we have forsaken these technical issues here. 18

1.0 0.8 0.6 ACF 0.40.2 -0.0 -0.2 0 20 Lag 40 60 Figure 5: The ACF plot for 16 first lags of the series Jt values given in Table 3 4.6 Bootstrapping the Surplus-Value Empirical Process In the sequel, let M denote a 10-th years surplus value of our company described in Section 4.4. In general, as described above M is a complicated function of a variety of parameters, some determin- istic and some stochastic in nature as well as the investment strategies. However, since in this work we are concerned with modeling cash flow as a function of underlying interest rates process for our purpose here we consider M to be a function of the time series {Jt } satisfying assumptions (A1) and (A2) above, assuming all other parameters to be fixed as described in Section 4.4. Under these conditions we may consider M to be a random variable with some unknown probability distribution H(x) = P rob(M ≤ x). Our goal is to obtain an estimate of H(x). In his paper KL fitted first the prescribed form of a probability distribution into the interest rates process and then used that model to approximate H(x) by generating interest rates scenarios via Monte-Carlo simulations. In contrast with his approach we shall proceed directly to the modeling of H(x) by bootstrapping the underlying empirical process. Let us consider a set of n = 167 values of Jt process obtained from annualized 30-year Treasury bond yield rates as given in Table 3. If we make an assumption that for the next 10 years the future values of increments of ln It will follow the same stationary process Jt , we may than consider calculating the values of M based on these available 167 data values. In fact, since the Jt ’s in Table 3 are based on the monthly rates, every sequence of 120 consecutive values of that series will give us a distinct value of M . This process results in n − 120 + 1 = 46 empirical values of M . However, by adopting the ”circular series” idea, i.e., by wrapping the 167 values from Table 3 around the circle and then evaluating M for all sequences of 120 consecutive Jt we obtain additional 119 empirical values of M and thus a total of n = 167 empirical values of M is available. The distribution of these values is presented in Figure 6. However, let us note that the obtained empirical M ’s are 19

Figure 6: 168 ordered values of the 10-th year surplus, each based on the 120 consecutive values of Jt series given in Table 3 using circular wrapping so that J1 = J167 , J2 = J168 , etc. not independent but rather form a time series which inherits the properties (A1) and (A2) of the b n (x) based original time series Jt . 4 A natural approximation for H(x) is now an empirical process H on the values of M obtained above, say, m1 , . . . mn (where in our case we would have n = 167), defined as X n b n (x) = 1 H I(mi ≤ x) (9) n i=1 where I(mi ≤ x) = 1 when mi ≤ x and 0 otherwise. By a standard result from the probability theory, under the assumptions that the mi ’s are realizations for a time series satisfying (A1) and (A2) we have that for any x Hb n (x) → H(x) as n → ∞. As we can see from the above considerations the empirical process H(·) b could be used to ap- proximate the distribution of random variable M and, in fact, a construction of such an estimator would typically be a first step in any non-parametric analysis of this type. In our case, however, the estimate based on the empirical values of the mi ’s suffers a drawback, namely that it uses only the available empirical values of surplus value, which in turn rely on very long realizations of the values of the time series Jt , while in fact explicitly using only very few of them. Indeed, for each mi we need to have 120 consecutive values of Jt , but in the actual calculation of surplus value we 4 This is not quite true, as there is perhaps a small contamination here due to the use of “wrapping” procedure for calculating empirical values of M . However, it can be shown that the effect of a circular wrapping quickly disappears as the number of available Jt0 s grows large. For the purpose of our analysis we have assumed the wrapping effect to be negligible. 20

use only 10 of them. This deficiency of Hb n (·) may be especially apparent for small and moderate values of m, i.e., when we deal only with a relatively short dependence structure of the underlying Jt process (and this seems to be the case, for instance, with our 30 year Treasury yields data in Table 3). Figure 7: The first four out of of 1000 generated distributions the bootstrapped 10-th year surplus value M ∗ . In view of the above, it seems reasonable to consider a bootstrapped version of the empirical process (9), say Hb n∗ , treated as a function of the underlying process Jt with some appropriately chosen block length l. Thus, our bootstrap estimator of H(x) ant any fixed point x (i.e., θ̂n∗ in our notation of Section 2.4, with H b n being θ̂n ) is X B b n∗ (x) = 1 H b ∗ (x) H (10) (b) B b=1 where H b b ∗ (x) is a value of H(x) calculated for the particular realization of the pseudo-values (b) J1 , . . . , Jn∗ generated via circular bootstrap procedure using k blocks of length l as described in ∗ Section 2.4, and B is the number of generated pseudo-sequences. 5 Let us note that the above estimate may be also viewed as a function of bootstrapped surplus values of M ∗ , say, m∗1 , . . . , m∗n obtained by evaluating the surplus value for pseudo-values of the series J1∗ , . . . , Jn∗ . For the purpose of the interval estimation of H(x) we consider bootstrap confidence intervals based on the method of percentiles, as described in Section 2.3 where now the bootstrap distribution G∗ denotes the distribution of H b ∗ (x) at any fixed point x, given the observed values of J1 , . . . , Jn . n 5 In order to achieve a desired accuracy of the bootstrap estimator for m-dependent data,the number of generated resamples is usually taken to be fairly large (typically, at least 1000). Of course, the number of all possible distinct resamples is limited by the number of available distinct series values (n) and a block length (l). 21

b Figure 8: Empirical surplus value process H(x) (thick line) and the four realizations of the boot- b n (x) calculated for the distributions given in Figure 7. The values strapped surplus value process H ∗ of x are given in millions. 4.7 Distribution of the 10–th Year Surplus. Comparison with the Para- metric Models In this section we shall discuss the results of our non-parametric cash-flow analysis as well as and give their brief comparison with the results obtained in KL. Using the cash flow model outlined above in Section 4.3 we analyze distribution of the 10-th year surplus using the data for average annualized yields to maturity for 30 years Treasury Bonds (as in KL, see Table 3). We have also presented the results of a similar analysis based on a larger set of data, namely on the average yield to maturity on long-term Treasury Bonds for the period 1953-1976. 4.7.1 Bootstrap Estimates of the Surplus Cumulative Distribution The bootstrap estimates H b n∗ (·) of the 10-th year surplus values distribution H(x) = P rob(M ≤ x) for a grid of x values ranging between (in millions) along with their 95% confidence intervals were calculated using a method of circular bootstrap as described in Section 2.4 with a circular block length of six months (i.e., in notation of Section 2.4 l = 6). Other values of l were also considered and the obtained results were more-less comparable for 1 < l < 6. In contrast, the results obtained for l ≥ 7 were quite different, however, they were deemed to be unreliable in view of a relatively small sample size (n=167), √ since the theoretical result on block bootstrap states that l should be of order smaller than n in order to guarantee the consistency of bootstrap estimates (cf. Shao and Tu, 1995 Chapter 9). The case l = 1 corresponds to an assumption of independence among the Jt ’s and was not considered here. 22

You can also read