The lifecontingencies Package: Performing Financial and Actuarial Mathematics Calculations in R

It is possible to model life contingency insurances with the lifecontingencies R package, which is capable of performing financial and actuarial mathematics calculations. Its functions permit one to determine both the expected value and the stochastic distribution of insured benefits. Therefore, life insurance coverage can be priced and portfolios risk-based capital requirements can be assessed. This paper briefly summarizes the theory regarding life contingencies that is based on financial mathematics and demographic concepts. Then, with the aid of applied examples, it shows how the lifecontingencies package can be a useful tool for executing routine, deterministic, or stochastic calculations for life-contingencies actuarial mathematics.


Introduction
This vignette is based on the package's paper published in JSS, Spedicato (2013). Apart of keeping track of package's updates, a major difference with respect to JSS publication is that parallel features are turned off to cope with CRAN packages' submission policies. As of March 2014, the lifecontingencies package (Spedicato 2013) appears as the first R package that deals with life contingent actuarial mathematics. The R statistical programming environment (Team 2012) has become the primary reference software for academics. Even in a business context, R is now considered a valid alternative to affirmed proprietary packages for statistics and data analysis, namely SAS (SAS Institute Inc. 2011), MATLAB (The MathWorks, Inc. 2011) and SPSS (IBM Corp 2012). Some packages for actuarial applications have already been developed within R. However, most of them mainly focus on non-life insurance. In fact, non-life insurance modeling involves more data analysis and applied statistical modeling than that of life insurance. Functions allowing one to fit loss distributions and perform credibility analysis are provided within the package actuar (Christophe Dutang, Vincent Goulet, and Mathieu Pigeon 2008). This package represents the computational side of the classical actuarial textbook on loss distribution (Klugman, Panjer, Willmot, and Venter 2009). The package ChainLadder (Gesmann and Zhang 2011) provides functions that are capable of estimating loss reserves for non-life insurance. Generalized linear models (GLMs), widely used in non-life insurance rate-making, can be fit by functions bundled within base R distributions. The generalized additive models for location, shape and scale (GAMLSS) and tweedie regression, which are both more advanced predictive models used by actuaries, are handled by specifically developed packages such as gamlss (Rigby and Stasinopoulos 2005;

Statistical and financial foundations of life contingencies
The actuarial pricing and reserving of life contingent insurances involves the calculation of statistics regarding occurrences and amounts of future cash flows. For example, the insurance pure premium (also known as benefit premium) can be regarded as the expected value of the prospective benefits cash flow distribution, valued at time zero for a given interest rate structure. The probabilities of the prospective benefits cash flow are based on the occurrence of the policyholder's life events (life contingencies). In addition, the theory of interest is used to find the present value of these amounts that will occur in the future. Therefore, life insurance actuarial mathematics is based on concepts derived from demography and theory of interest.
A life table (also called a mortality table or actuarial table) is a table that shows how mortality affects subjects of a cohort across different ages. For each age x, it reports the number of l x individuals living at the beginning of age x. It is a sequence of l 0 , l 1 , . . . , l ω , where l ω , the terminal age, represents the latest age that a subject of the cohort can survive until. Life tables are typically distinguished according to gender, year of birth and nationality, with different categories being used depending on the type of life contingency (i.e., assurance or annuity). As another example, businesses may have different customers with different underlying mortalities, so they will be in need of different life tables.
From a statistical point of view, a life table allows one to deduce the probability distribution of the future lifetime for a policyholder aged x. In particular, a life table also permits one to derive two key probability distributions:T x , the complete future lifetime for a policyholder aged x, and its curtate form,K x , the number of future years completed before death. Therefore, many demographic statistics can be derived from a life table, of which a non-exhaustive list follows: • t p x = l x+t lx , the probability that a policyholder alive at age x will reach age x + t. • t q x , the complementary probability of t p x .
• t d x , the number of deaths between age x and x + t.
• t L x = t 0 l x+y dy, the expected number of years lived by the cohort between ages x and x + t.
• t m x = tdx tLx , the central mortality rate between ages x and x + t.
• e x , the curtate expectation of life for a subject aged x, e The Keyfitz textbook (Keyfitz and Caswell 2005) provides an exhaustive coverage of life table theory and practice. Life tables are usually published by institutions that have access to large amounts of reliable historical data, like social security bureaus. It is common practice for actuaries to start from these life tables and adapt them to the insurer's portfolio actual experience.
Classical financial mathematics deals with monetary amounts that can be available at different times. The present value of a series of cash flows, expressed by Equation 1, is probably the most important concept. The present value (PV) can be considered as the value -in current money -of a series of financial cash flows, CF t , that are available at different periods of time.
The interest rate, i, represents the measure of the price of money available in future times. Like the interest rate, the time value of money can be expressed by the discount rate d = i 1+i . This paper will use the i symbol to express effective compound interest, with money invested once per period. In the case where money is invested more frequently, say m times per period, each fractional period is named the interest conversion period. During each interest conversion period, the real interest rate i (m) m is earned, where the i (m) expression defines the convertible (also known as "nominal") rate of interest payable m times per period.
Equation 2 combines the various notations for interest and discount rates, both on an effective and convertible basis, to express what an amount of $1 invested now will grow to at time t.
All financial mathematics functions (like annuities, a n , or accumulated values, s n ) can be rewritten as particular expressions of Equation 1, as shown in classical actuarial mathematics textbooks.
Actuaries use the probabilities inherent in life tables to price life contingent insurances. In fact, life contingencies are stochastic variables themselves. A life contingent insurance can be represented by a sequence of one or more payments whose occurrence, timing, and consequent present value are not certain. In fact, their eventual occurrence and timing depend on events in the life of the policyholder; for this reason, they are named "life contingencies". Since the actuary focuses on the present value of such uncertain payments, future payments of life contingent insurances need to be discounted using interest rates that may also be considered stochastic. The lifecontingencies package provides functions that model most of the standard life contingent random variables,Z, and in particular their expected value, the actuarial present value (APV). Of all the statistics used by actuaries, the APV is certainly the most important. In fact, it represents the average cost of the benefits guaranteed to the policyholder by the insurer. In a non-life insurance context, it would be named pure premium. The policyholder pays out the gross premium, G, which is a sum of benefit premiums, loading for expense, profits and taxes. Life contingencies can be either continuous or discrete, as the cited actuarial mathematics textbooks explain. Directly, the lifecontingencies package can only model discrete life contingencies with a non-stochastic interest rate. Nevertheless, most continuous time life contingent insurances can be easily derived from their discrete form under broad assumptions that can be found in the cited textbooks.
A few examples of life contingent insurances follow: 1. An n-year term life insurance provides a payment, if the insured dies within n years from issue. If the payment is performed at the end of the year of death,Z can be written asZ Its APV expression is A 1 x:n . 2. A life annuity consists of a sequence of benefits paid out as long as the insured life survives. In particular, a temporary life annuity due pays a benefit at the beginning of each period as long as the annuitant aged x survives, for up to a total of n years, or n payments.Z can be written asZ = ä K+1♣ ,K x < n, a n♣ ,K x ≥ n.
Its APV expression isä x:n .
3. An n-year pure endowment insurance grants a benefit payable at the end of n years if the insured survives at least n years from issue.Z can be written as Its APV expression is A 1 x:n (or n E x ). 4. An n-year endowment insurance will pay a benefit at the year of death or at the end of the n-th year, whichever occurs earlier.Z can be written as Its APV expression is A x:n .
Interested readers could see the cited references for formulas regarding other life contingent insurances like (DA) 1 x:n , the decreasing term life insurance, or (IA) 1 x:n , the increasing term life insurance. There are also common variations of payments that form arrangements like deferment or fractional payments. Similarly, it is possible to define insurances and annuities depending on the survival status of two or more lives. For example, A xy and a xy represent, respectively, the APV symbols for the two lives joint-live insurance and the two lives lastsurvivor annuity immediate.
The lifecontingencies package provides functions that allow an actuary to perform classical financial and actuarial mathematics calculations. In addition to standard deterministic modeling, a peculiar feature of lifecontingencies is that it allows one to generate variates from the stochastic distribution of the present value of future benefits,Z, for most life contingent insurances. This feature allows for a deeper assessment of the insurance liabilities variability.

The structure of the package
The package lifecontingencies contains classes and methods that handle life tables and actuarial tables in a convenient manner.
The package is loaded within the R command line interace as follows:

R> library("lifecontingencies")
Two main S4 classes have been defined within the lifecontingencies package: the 'lifetable' class and the 'actuarialtable' class. The 'lifetable' class is defined as follows:
Finally, the lifecontingencies package depends on the methods package, which defines its classes, and the parallel package, which is used to speed up computations. As detailed in Section 5, implementation of C or C++ code snippets is expected to shorten computational times in the future versions of the package.

Code and examples
This section is structured as follows: Section 4.1 shows classical financial mathematics examples, Section 4.2 deals with life tables and actuarial table management, Section 4.3 shows classical actuarial mathematics examples, and Section 4.4 presents functions in the lifecontingencies package that perform simulation analysis.

Classical financial mathematics example
The lifecontingencies package provides functions that perform classical financial mathematics Function Purpose presentValue present value of a series of cash flows annuity present value of an annuity-certain, a n accumulatedValue future value of a series of cash flows, s n increasingAnnuity present value of an increasing annuity -certain, IA n decreasingAnnuity present value of a decreasing annuity -certain, DA n convertible2Effective conversion from convertible to effective interest (discount) rates effective2Convertible convertible2Effective inverse intensity2Interest conversion from force of interest to the interest rate interest2Intensity intensity2Interest inverse duration dollar / Macaulay duration of a series of cash flows convexity convexity of a series of cash flows  Table 2.
Some of these implement closed form formulas and their inverses; this is also shown in financial mathematics textbooks. A broader discussion, however, shall be dedicated to the presentValue function. In fact, the presentValue function is internally called by most other financial and actuarial functions within the lifecontingencies package. This function calculates present values or APVs by computing Equation 3.
The terms in Equation 3 are the cash flows vector, c i , the corresponding discount factors vector, v t i , and the occurrence probabilities vector, p i , respectively. Many lifecontingencies package functions, like axn or annuity, work by first defining the pattern vectors of cash flows, interest rate and probabilities (in case of actuarial functions), which are then passed as arguments to the presentValue function.
Examples that follow show how to handle interest and discount rates with different compounding frequencies, how to perform present value, annuity, and future value analysis, and how to deal with loan amortization and bond pricing.

Interest rate functions
Interest rates represent the time-value of money. Different types of rates can be found in literature. As a remark, Equation 4 displays the relationship between the effective interest rate, the convertible interest rate, the discount factor, the force of interest, the effective discount rate and the convertible discount rate. ( The functions interest2Discount, discount2Interest, convertible2Effective, effective2Convertible, interest2Intensity, intensity2Interest have all been based on Equation 4; their inverse formulas are implied therein. Throughout the paper, an effective interest rate is used unless otherwise stated.
As examples, functions interest2Discount and discount2Interest represent a convenient way to switch from interest to discount rates and vice versa.
[1] 0.03 The function convertible2Effective allows one to find the effective interest rate implied in a consumer -credit loan that offers a 10% convertible (nominal) interest rate with quarterly compounding.

Analysis of present value and internal rate of return
Performing a project appraisal means evaluating the net present value (NPV) of all projected cash flows. The code below shows an example of NPV analysis.

Annuities and future values
An annuity (certain) is a sequence of payments with a specified amount that is present valued. If it is valued at the end of the term of payment is is called future value (or accumulated value). The code below shows examples of annuity, a n♣ , and accumulated value, s n♣ , evaluations. The PV of an annuity immediate $100 payable at the end of each year for the next 5 years at an interest rate of 3% is: R> 100 * annuity(i=0.03, n=5) [1] 457.9707 while the corresponding future value is: R> 100 * accumulatedValue(i=0.03, n=5) [1] 530.9136 Annuities and future values payable k-thly (where fractional payments of 1 k are received for each k-th of a period) can be evaluated properly by setting the functions' parameters as follows: R> ann1 <-annuity(i=0.03, n=5, k=1, type="immediate") R> ann2 <-annuity(i=0.03, n=5, k=12, type="immediate") R> c(ann1,ann2) [1] 4.579707 4.642342 increasingAnnuity and decreasingAnnuity functions handle increasing and decreasing annuities, whose symbols are IA x , DA x respectively. Assuming a ten-year term and a 3% interest rate, examples of increasing and decreasing annuities follow.

Loan amortization
As this section exemplifies, lifecontingencies financial mathematics functions allow one to define the repayment schedule of any loan arrangement. Let C denote the loaned capital (principal). Assuming an interest rate i, the amount due to the lender at each installment is R = C a n| . Therefore, the R amount repays I t = C t−1 * i as interest, while C t = R − I t is the amount of the loan still outstanding after installment t has been paid. The periodic installment of loan repayment , R, is calculated as follows: R> capital <-100000 R> interest <-0.05 R> payments_per_year <-2 R> rate_per_period <-(1+interest)^(1/payments_per_year)-1 R> years <-30 R> R <-1/payments_per_year * + capital/annuity(i=interest, n=years, + k=payments_per_year) R> R [1] 3212.9 Then the balance due at end of the period (EoP) is calculated as follows: Figure 1 shows the EoP balance due for a 30 year loan, assuming a 5% interest rate on a principal of $ 100,000.

Bond pricing
Bond pricing represents another application of present value. A standard bond with face value C and term length T consists of equal coupons c paid at regular intervals. The final payment at time T is C T + c. Equation 5 expresses the present value of a bond with n remaining coupons.
Perpetuities are financial contracts that offer an indefinite sequence of payments either at the end (perpetuity-immediate) or at the beginning of each period (perpetuity-due). Loan amortization year EoP balance due As displayed below, the bond and perpetuity functions defined above can be used to price any bond, given face value, coupon rate, and term.

Duration and ALM
As defined within the package, duration and convexity formulas are reported in Equation 6 and Equation 7 respectively. Their typical application lies within porfolios' asset -liability management (ALM). The interested reader can find details in Chris Ruckman and Joe Francis (2006) and Broverman (2008) textbooks. However, the following example shows how the Macaulay duration (ex1), modified duration (ex2), and convexity (ex3) of any series of cash flows can be calculated by lifecontingencies package functions.
2. A perpetuity-immediate. As a remark, the formulas for the PV, duration and convexity of a perpetuity immediate are P V pp = 1 y , D pp = 1+y y , C pp = 2 y 2 respectively, if the yield rate is y.
Assume the issuing company wants to hedge its liability with an investment portfolio that will not be affected adverserly by changes in the investment yield. In order to solve the ALM problem, the composition of assets within the portfolio shall be chosen accordingly. Moreover, assume that the current market yield rate is 4%. The following lines of code figure out some parameters that are used within the example.
2. Setting the interest rate sensitivity (i.e., the duration) of assets to be equal to the interest rate sensitivity of liabilities. This is done by solving the system of equations shown in Equation 8. The parameters w i and D i stand for asset hedging weights and duration values respectively.
3. Setting the convexity of assets to be greater than the convexity of liabilities. In other words, this means verifying that asset decline (growth) will be slower (faster) than liability decline in case of a change in the interest rate.
The following lines of code calculate the asset weights vector by linear algebra functions bundled in R base.
R> a <-matrix(c(durBond, durPpty,1,1), nrow=2, + byrow=TRUE) R> b <-as.vector(c(7,1)) R> weights <-solve(a,b) R> weights [1] 0.8924163 0.1075837 Vector weights displays the portfolio composition in terms of bonds and perpetuities, respectively. Therefore, the number of bonds and perpetuities that can be purchased is determined by: It can be verified that the convexity of assets is greater than the convexity of liabilities.
R> convAsset <-weights[1] * convBond + weights[2] * covnPpty R> convAsset>convLiab The portfolio is immunized from yield rate variations because the present value of assets will be greater than the present value of the liabilities if the interest rate suddently drops to 3% just after hedging the asset purchase. The same occurs in case of an upward shift in the interest rate toward 5%.

Analysis of life tables and actuarial tables
Function Purpose dxt deaths between age x and x + t, t d x . pxt survival probability between age x and x + t, t p x . pxyzt survival probability for two (or more) lives, t p xy . qxt death probability between age x and x + t, t q x . qxyzt death probability for two (or more) lives, t q xy . Txt number of person-years lived after exact age convert death probabilities into mortality rate. mx2qx convert mortality rate into death probabilities. exn expected lifetime between age x and age x + n, n e x . rLife sample from the time until death distribution underlying a life table. rLifexyz sample from the time until death distribution underlying two or more lives. exyz n-year curtate lifetime of the joint-life status. probs2lifetable life table l x from raw one -year survival / death probabilities. lifetable and actuarialtable classes are designed to handle demographic and actuarial mathematics calculations. An actuarialtable class inherits from lifetable class; it adds one more slot for the rate of interest. Both classes have been designed using the S4 R classes framework. Table 3 lists the functions that have been developed for performing demographic analysis within lifecontingencies package. This section briefly exemplifies these functions.

Creating lifetable and actuarialtable objects
Life table objects can be created by using either raw R commands or existing data.frame objects. However, three components are needed to build a lifetable class object: 1. The years sequence, which is an integer sequence 0, 1, . . . , ω. It shall start from zero and end at ω, the terminal age (the age x for which p x = 0).
2. The l x vector, which is the number of subjects living at the beginning of age x; in other words, the number of subjects at risk of dying between year x and x + 1.
3. The name of the life table.
There are three main approaches for creating a lifetable object: 1. Directly from the x and l x vector.
2. By importing x and l x from an existing data.frame object.
3. From using raw survival probabilities.

R> head(exampleLt)
x lx 1 0 1000 2 1 950 3 2 850 4 3 700 5 4 680 6 5 600 Still, the easiest way to create a lifetable object is to start from a suitable existing data.frame. This will probably be the most practical approach for practicing actuaries. Some life or mortality rate tables have been bundled within the lifecontingencies package, as  The following example shows how US Social Security life tables are loaded from the existing demoUsa data set bundled in the lifecontingencies package.

R> getOmega(exampleAct)
[1] 9 Method print behaves differently for lifetable and actuarialtable objects. In fact, when the print method is applied on a lifetable object, it tabulates both the one year survival probability and the complete expected remaining life until death. Conversely, when the print method is applied on a lifetable object, classical commutation functions (D x , N x , C x , M x , R x ), discussed later, are printed out.

R> print(exampleLt)
Life It is possible to convert the actuarialtable object into a data.frame object, as shown below.
R> exampleActDf <-as(exampleAct, "data.frame")  A recent lifecontingencies package enhancement allows to export a life table as non -homogeneous discrete Markov chain by means of markovchainList S4 class object as defined in markovchain (Spedicato, Giorgio Alfredo 2015) R package.
R> data(soa08) R> require(markovchain) R> soa08Mc<-as(soa08,"markovchainList") Finally, a plot method can be applied to eitherlifetable or actuarialtable objects. The underlying survival function (which is the plot of x vs l x ) is displayed in both cases. Figure 2 shows the plot methods applied on a Society of Actuaries (SOA) actuarial table at 6% interest, which comes bundled within the lifecontingencies package as soa08Act object.

Basic demographic analysis
Basic demography calculations can be performed on valid lifetable or actuariatable objects. The functions discussed in this section calculate proper ratios or sums on l x or d x values using functions that access to the lifetable object slots. This is all done in accordance with the demographic formula definitions.
R> data("soa08Act") R> pxtLin <-pxt(soa08Act,80,0.5,"linear") R> pxtCnst <-pxt(soa08Act,80,0.5,"constant force") R> pxtHyph <-pxt(soa08Act,80,0.5,"hyperbolic") R> c(pxtLin,pxtCnst,pxtHyph) [1] 0.9598496 0.9590095 0.9581701 Calculations for survival probabilities on two (or more) lives can also be performed. As a remark, two different life statuses are defined within the analysis of multiple lives survival: "joint" survival status and "last" survival status. The "joint" survival status exists while all the members of the pool are alive, while the "last" survival status exists until the last member of the pool dies. All calculations assume that the multiple lives are independent. Equation 9 expresses the remaining future lifetime on a couple xy under the joint and last survival status respectively.
T xy = min (T x , T y ) The following code shows how the joint survival probability, last survival probability, and expected joint lifetime can be evaluated using functions in lifecontingencies.

Life insurance examples
The evaluation of the APV has traditionally followed one of three approaches: the use of commutation tables, the current payment technique, or the expected value method. Commutation tables extend the life table by tabulating special functions of age and rate of interest, as Anderson (1999) further considers. Ratios of commutation table functions allow an actuary to evaluate APV for standard insurances. However, commutation table usage has become less prominent in the computer era. In fact, these tables are not flexible enough and their usage is computationally inefficient. Therefore, the lifecontingencies package does not use the commutation table approach to evaluate APVs.
The current payment technique calculates the APV of a life contingency insurance,Z, as the scalar product of three vectors:Z = ⟨⟨c •v⟩ •p⟩; this uses the vector of all possible uncertain cash flows,c, the vector of discount factors,v, and the vector of cash flow probabilities,p. The lifecontingencies package implements the current payment technique by using actuarial functions listed in Table 5 to evaluate APVs. Finally, the expected value approach modelsZ as the scalar product of two vectors:Z = p k •x .pk is P r K = k , the probability that the future curtate lifetime will be exactly k years, wherex is the present value of benefits due under the policy term ifK = k. rLifeContingencies and rLifeContingenciesXyz implement the expected value approach to generateZ variates.
Consider an n year annuity due. Its APV,ä x:n , using the commutation table approach is reported in Equation 10, while Equation 11 reports the same APV using the current payment technique. Finally, Equation 12 calculates the APV using the expected value approach.
In order to understand how lifepackage implements the current payment technique in its actuarial function, it is worthwhile to look closer at the core of the axn function. This function takes the following parameters as inputs: n, the term of the annuity; k the fractional payment frequency; x the annuitant age; m, the deferring period. Then, it defines: 4. Finally, the three vectors are passed as input parameters to the presentValue function as the following code shows: presentValue(cashFlows=payments, timeIds=times, interestRates = interest, probabilities=probs) In the examples that follow, the SOA illustrative actuarial table is used in the calculation of premiums and reserves of life contingencies.
The first example values a 40-year insurance on a policyholder aged 25, with benefits payable at the end of the month of death. Equation 13 would determine the benefit premium using the commutation table approach.
The following lines of code compute the benefit premium using UComm, the commutation technique, and UCpt, the current payment technique. If, while the policyholder is alive, the premium is paid in ten equal installments at the beginning of each year instead of a lump sum, then the yearly premium, P , would be determined as follows: R> P <-UCpt/axn(actuarialtable=soa08Act,x=25,n=10) R> P [1] 0.006351046 The lifecontingencies package allows one to evaluate APVs of endowment insurances; this can be calculated for increasing and decreasing life insurances as well. The lines of code that follow will prove the actuarial equivalence expressed by Equation 14 in a computational context.

Life annuity examples
Life contingent annuities form sequences of payments whose occurrence and duration depend the policyholder's future lifetime. The few examples that follow demonstrate how the lifecontingencies package can directly compute the APV for typical life contingencies using either bundled functions or classical commutation tables.

Benefit reserves examples
The (prospective) benefit reserve consists in the difference between the APV of obligated future benefit payments due by the insurer and the APV of the projected inflows due by the policyholder. It represents the outstanding net insurer's obligation arising from the underwritten insurance policy. An example will clarify this concept. The code below evaluates the benefit reserve for a 25 year old 40 -year life insurance of $ 100,000, with benefits payable at the end of the year of death, assuming a level benefit premium payable at the beginning of each year. The benefit premium and reserve equations for this life contingent insurance are displayed by Equation 17. R> P=100000 * Axn(soa08Act,x=25,n=40)/axn(soa08Act,x=25,n=40) R> reserveFun = function(t) return(100000*Axn(soa08Act,x=25+t,n=40-t)-P* + axn(soa08Act,x=25+t,n=40-t)) R> for(t in 0:40) {if(t%%5==0) cat("At time ",t, + " benefit reserve is ", + reserveFun(t),"\n")} Another reserve calculation example shows the benefit reserve for a deferred annuity-due on a policyholder aged 25 when the annuity is deferred until age 65. The code below shows the reserve calculation while Figure 3 plots the outstanding reserve at the end of each contract year.

Expense considerations
The premium paid by the policyholder usually contains an allowance for expenses and profit loading. Expenses cover policy servicing and the producers' commissions. The insurers' profit load is explicitly taken into account in the benefit premium as a flat amount, or, sometimes, as a percentage of the final premium. In other cases implicit profit loading is generated by using demographic and financial assumptions more prudentially than would be necessary. The equivalence principle can be extended to both the gross premium, G, and the expense augmented reserve, t V E , when expense allowances are taken into account by using Equation 18.
The following example shows how an expense loaded premium G is calculated for a $ 100,000 whole life insurance on a 35 year old policyholder. 10% of premium expense per year, 25 policy expenses per year, and annual maintenance expense of 2.5 per 1,000 units of capital are assumed.

Insurances and annuities on two lives
The package provides functions designed to evaluate life insurances and annuities on two lives. The following example checks the actuarial mathematics identity on joint and last survival status annuities expressed by Equation 19.

Stochastic analysis
This last section illustrates some stochastic analysis that can be performed by the lifecontingencies package, in both demographic ( Section 4.4.1 ) and life insurance contexts ( Section 4.4.2 ).

Demographic examples
The age-until-death, both in the continuous,T x , or curtate form,K x , is a stochastic variable whose distribution is intrinsic in the deaths within a life table. Therefore, a dedicated function, rLife, has been designed within the lifecontingencies package to draw a sample from either K x orT x . Drawing from K x is quite simple: the distribution of the curtate future lifetime is defined as Pr K , and it is passed as a prob parameter to base R sample function. For example, the code below shows how the rLife function can be used to draw a sample of size five from the curtate future lifetime of a policyholder aged 45 as implied by the SOA life table.
R> rLife(n = 5, object = soa08Act, x = 45, type = "Kx") [1] 40 18 29 12 51 rLifexyz represents the multiple lives extension of the rLife function. It returns a matrix of sampled expected future lifetimes of J policyholders given a list of J lifetables. The simulation approach is useful in evaluating demographic quantities when the analytical approach is not feasible. One example could be the expected years of widowhood, which Equation 20 defines. T x andT y in Equation 20 stand for complete future lifetimes for the husband and the wife respectively.
The following code shows how this function could be used to evaluate the expected years of widowhood for the wife of a couple. The example makes use of the Italian projected life tables ips55M and ips55F, whose derivation was shown in Section 4.2.

Distribution of widowance yars
[1] 7.43 Finally, Figure 4 shows the distribution of widowance years as determined in the previous example.

Actuarial mathematics examples
The present value of the future benefits cash flows distribution,Z, is a random variable. It is a function of both the interest rate and the indicator variables which are determined by the life status of the insured. Both of these quantities can be deemed stochastic. However, interest rates are considered deterministic within the framework of the current version of lifecontingencies package. The generation of n-size variates fromZ is performed by the following algorithm: 1. Define a function, P V that returns the present value of the life contingent insurance benefits, given the age at death of the policyholder, as T 0 , P V (T 0 ). Within the lifecontingencies package, present value functions have been defined for the most important life contingencies. Such functions are not visibly exported in the package namespace.
3. Use T 0 variates as inputs for P V (T 0 ) to get variates fromZ.
The code below shows the internal function .faxn, which returns the present value of a life contingent insurance. .faxn is internally called by the rLifeContingencies function, as discussed below. T, y, n, i, m, k represent the age at death, the attained age, the term of the annuity, the interest rate, the deferring period, and the fractional payment frequency respectively.
.faxn<-function(T,y,n, i, m, k=1) { out <-numeric(1) K <-T-y if(K<m) { out <-0 } else { times <-seq(from=m, to=min(m+n-1/k,K),by=1/k) out <-presentValue(cashFlows=rep(1/k, length(times)), timeIds=times, interestRates=i) } return(out) } Life contingency insurance functions return the APV, E Z , as a default value. The functions in Table 5 compute APVs by using the current payment technique. Another possible approach for evalutating APVs, even if computationally inefficient, could be to draw a sample from the underlyingZ distribution and compute its sample mean. Table 5 returns a sample of size one if the type parameter default value, "EV" ( expected value), is overridden by the string "ST" ( stochastic).

Every function in
However, when samples of greater size are required, the most straightforward approach is the rLifeContingencies function. The code below shows how to generateZ variates from either term life insurances, increasing term insurances, temporary annuities, or endowment insurances respectively. For each example, the lack of bias is verified by comparing the mean of the sample with the theoretical APV using a classical t -test. All examples refer to an individual aged 20 with a 40 year insurance. Figure 5 shows the resultingZ distributions.
The full distribution of a life contingent insurance variableZ, can be used to compute premiums using the percentile premium principle. Under this approach, the premium is set to ensure that the insurer will suffer financial loss with a sufficiently low probability (made explicit by the percentile). An example will clarify the concept. For a 40 -year insurance on a single policyholder aged 25, the actuarial present value of benefits, i.e., the expected value of discounted future benefits, would be R> APV <-Axn(actuarialtable = soa08Act, x=25, n=40) R> APV [1] 0.04797088 while the benefit premium at the 90th percentile, that is, the premium that would make the insurer incurr an underwriting loss with 10% probability, would be R> samples <-rLifeContingencies(n=numSim, lifecontingency = "Axn", + object= soa08Act, x=25,t=40,parallel=FALSE) R> pct90Pr <-as.numeric(quantile(samples,.90)) R> pct90Pr [1] 0.1234772 Finally, if N = 1000 similar policyholders were insured, the Law of Large Numbers would lead to a strong reduction in the premium charged on each policyholder, as computed below.
R> pct90Pr2 <-qnorm(p=0.90,mean=APV, sd=sd(samples)/sqrt(1000)) R> pct90Pr2 [1] 0.05307155 The final example in this paper shows how the stochastic functions bundled in the lifecontingencies package can be used to make an actuarial appraisal of embedded benefits. Suppose a corporation grants its 100 employees life insurance benefits equal to their annual salary and payable at the month of death. Moreover, suppose that: 1. The expected value and the standard deviation of the salary are $ 50,000 and $ 15,000 respectively and the salaries follow a log-normal distribution.
2. The distribution of the employees' age is uniform on [25,65]. Assume 65 is the retirement age.
3. The SOA illustrative table represents an unbiased description of the population mortality.
4. Assume no lapse to hold.

The policy length is annual.
We evaluated the best estimate, or the fair value of the insured benefits according to both IFRS accounting standards and risk margin measure. In this example, the risk margin measure, which is the difference between the 75th percentile and the best estimate, will be used. IFRS standards (Post, Grandl, Schmidl, and Dorfman 2007) define the fair value of an insurance liability as the sum of its best estimate and its risk margin.
In the initial part of the example, we set up the parameter of the model and configured the parallel computation facility available with the package parallel.

Multiple Decrement Models within lifecontingencies Package
Until now no R package provides a good tool to manage multiple decrement tables, even if Deshmukh (2012) provides an R based focus on multiple decrement tables with applications in R. The topic is deeply related to multistate analysis of life histories on which Willekens (2014) provide a very good introduction. The lifecontingencies package's mdt class that has been specifically engineered to manage multiple decrements models with R. Applied examples will follows.
Following notation in Finan (2014), we provide definitions of the key quantities that allow to understand the main concepts regarding Multiple Decrement (MD) theory. Be l x survivors to age x that will, at future ages, be fully depleted by m causes of decre- represents the expected number of lives exiting from the population between ages x and x + 1 due to decrement j. Therefore it follows that n d The probability that a life x will leave the group within one year as a result of decrement j x and that t q x .

The mdt class
Examples in this paper are worked on slides provided in Valdez (2011).
We create a mdt class object. We can use the first example found on (?, p. 4 The mdt class is an S4 class object (Chambers 2008) comprised by a character slot name and a data.frame slot table that is composed by following columns: 1. x: the age, from 0 to ω.
2. lx: the subject living (at risk) at the beginning of age.
3. one or more colums for different causes of decrements.
Values within table item represents absolute number of subjects at risk at the beginning of age x and dying for cause j during period x -x + 1.
Within the various methods defined within the mdt class, setValidity performs consistency checks to properly create the mdt object. In particular, it verifies whether: 1. x and lx exist and that they are consistent. x should start from 0 and flows by increments of one. The first lx value should be equal to the sum of all decrements and that 2. If the decrements (or x and lx) have been provided only for partial ages, the table is completed below (from 0 to l x−1 ) assuming a decrement rate of 0.01 for the first cause of death.
As shown, when the table is sanitized the operations performed are reported on logs.
An internal function, .tableSanitizer tries to fix the limitations on the input table in order it to meet the class definition requirements. Similarly, it is possible to export a mdt to a data.frame or to a markovchainList object (from markovchain package).

R> print(valdezMdt)
R> valdezDf<-as(valdezMdt,"data.frame") R> require(markovchain) R> valdezMarkovChainList<-as(valdezMdt,"markovchainList") Two specific methods have been defined for mdt class objects: getOmega, that returns the maximum attainable age (similar to the one of lifetable class), and getDecrements, that returns the decrements (by means of the names within table slot different from x and lx).

R> summary(valdezMdt)
This is Multiple Decrements

Decrement probabilities calculation
The lifecontingencies package makes easy to compute d x as well as n d (τ ) x quantities thanks to dxt function.

Associated Single Decrement Table calculation
For each force of decrement µ j (x + t) the Associated Single Decrement Table (ASDT) is a decrement models that assumes survivoship depends only from j. Within ASDT the equations shown in Equation 21 are true : Assuming uniform distribution of decrements (UDD), that is t q (j) x , t ∈ [0, 1], then Equation 22 is valid: This has been implemented in the qxt.prime.fromMdt function: R> qxt.prime.fromMdt(object = valdezMdt,x=53, decrement="accidents") [1] 0.0003504636 If UDD holds also for each associated single decrement Equation 24 is also valid: The Equation ?

Actuarial Applications
The package now offers limited capabilities to fit multiple decrement insurances, e.g. (A 1 x:n ) (1) The example in (Finan 2014, p. 674), cites: A 3-year term issued to (16) pays 20,000 at the end of year of death if death results from an accident. The mdt table is below created.

Completed the table at top, all decrements on first cause
The value of (A 1 16:3 ) (a) is below calculated R> Axn.mdt(object=myMdt,x=16,i=.1,decrement="da") [1] 0.1363636 Another example has been inspired by data found in De Angelis, Paolo and Di Falco, L. (2016). We use life tables from the quoted book and the following functions (still based on single decrement tables) to compute the APV of (costant or variable) annuity with multiple decrement. The first type of annuity pays benefits are payable while the annuitant is in current state and cease upon transition to another state.

Advantages, limitations, and future perspectives
The lifecontingencies package allows actuaries to perform demographic, financial and actuarial mathematics calculations within R; in particular, life contingent insurance contracts can be priced and reserved. In addition, a peculiar feature of lifecontingencies is its ability to generate variates from the future life time and the underlying stochastic distributions of life contingent insurances.
One of the most significant limitations of the most recent package release is that few functions are available to model multiple decrements tables. In addition, continuous-time life contingent models are currently not handled explicitly.
We expect to remove such limitations in the future. In addition, we expect to provide coerce methods for packages that specialize in demographic analysis, like demography and LifeTables. Furthermore, we wish to allow easier sharing of analyses with interest rate modeling packages like termstrc.
Finally, code optimization and improvement is carried out continuously. The extension of parallel computation features, memory usage profiling, and the use of C or C++ code fragments in select parts of the code have begun (for Rcpp package, Eddelbuettel (2013b) usage) or being planned for the near future.

Accuracy
The accuracy of the calculations has been verified by checking numerical examples reported in Bowers et al. (1997) and in the lecture notes of "Actuarial Mathematics" taken by the author years ago at The Catholic University of Milan (Mazzoleni 2000). Such test have been implemented with unit root testing in the package testthat, Wickham (2011). The numerical results are identical to those reported in the cited references for most functions, with the exception of fractional payment annuities; for these, the answer is only reported up to the fifth decimal. The reason for such inaccuracy is that the package calculates the APV using the direct sum of fractional survival probabilities, while the results reported in the cited references are obtained by closed formulas.
Finally, it is worth noting that the package and functions herein are provided as is without any guarantee regarding the accuracy of calculations. The author disclaims any liability arising from potential losses due to direct or indirect use of this package. Users are invited to notice the author any potential bug from the github package site https://github.com/ spedygiorgio/lifecontingencies.