AnnuityRIR: an R-package to approximate the value of an annuity according to the non-central moments of the capitalization factor

The aim of this paper is to design the package of the R statistical software called “Annuity Random Interest Rate”, referred hereinafter as AnnuityRIR, in order to calculate the value of an n-annuity with payments of one unit each when the interest rate is random. To do this, we have employed different approaches; the two main methodologies treated in this study consider that all non-central moments of the capitalization factor are known, or contrarily some of them are unknown. Consequently, five different approaches have been developed and the practical application of the proposed methods is reflected in this paper by pricing an annuity with a random risk-free interest rate during the last ten years. The version is available from CRAN: https://cran.r-project.org/web/packages/AnnuityRIR/index.html.


Introduction
An annuity is a sequence of n payments separated by the same interval of time (Gerber 1997).The appraisal of an investment, or the development of loan repayments schemes, pensions, bonds, mortgages, house rents, insurances and other derivatives contracts are just a few examples of these periodic payments known as annuities (Date et al. 2007).In the assessment of the present or the final value of an annuity, traditional actuarial theory has assumed that the interest rate is fixed for all the involved years (Bowers et al. 1997).In effect, classic literature on pricing annuities concentrates on using a deterministic discount factor (Kellison 1991).However, this approach is not valid for those operations embedded in uncertain environments, given that the interest rate is neither constant nor known (Burnecki et al. 2003).Therefore, it seems reasonable to let interest rates used in future valuations , vary as a random variable over time (Kellison 1991).As indicated by Duffee (2002), "the standard class of affine models produces poor forecasts of future Treasury yields, better forecasts are generated by assuming that yields follow random walks".
Nevertheless, in this paper, we will not derive the interest rate as a stochastic process over time but we will use the distribution function of the random variable obtained starting from the current values exhibited by the interest rate.In this way, the employment of interest rates as a part of the stochastic discount factors allows a more accurate pricing of any operation (Date et al. 2007;Dufresne 2007).Our proposal is that the calculation of the present and the final value of an annuity with random interest rate must be based on the real data corresponding to the EURIBOR at the present moment.Starting from these empirical , discrete values, we can determine easily the non-central moments of all (positive and negative) orders.But, if the data fit a well-known statistical distribution, we can take advantage of the accuracy of the distribution and then to determine the moments by using all the powerful information contained in these distributions.
In general, the assessment of operations in which some parameters are random variables requires a treatment through the formulation of potential scenarios, which are subsequently reduced to one by statistical treatment (Cruz Rambaud and Valls Martínez 2002).In the context of annuities assessment, the interest rate has a great relevance because even small changes in the rate may cause major changes in the total annuity value.Thus, the determination of the value of the interest rate should be carried out as accurately as possible.In this respect, many experts have discussed on the calculation of annuities with random rates of interest.In particular, Liu et al. (2011) study the present value of increasing, decreasing and paused rainbow immediate annuities with random interest rate.Specifically, we highlight Zaks' (2001) contribution, complemented by Burnecki et al. (2003) who modifies some Zaks' formulas.Both studies analyze deeply the formulation of annuities with random interest rates.As indicated by Baker (2001), in these studies, the expected value of a payment of $1 at the end of year k is represented as: being i k the rate of interest for year k = 1, 2, . . ., n.So, , it can be deduced that: being where the variance of the amount accumulated during the year k may be determined as var(1 (2) Zaks ( 2001) and Burnecki et al. (2003) focus on the mean and the variance of the final value of an n-year unitary annuity-due, (C n ), assuming that the interest rates i 1 , i 2 , . . ., i n are independent, and have the same mean, 1 + j, and variance, s 2 .In this case, the following identity holds (Zaks 2001, Theorem 3.1).
Unfortunately, the assumption of independent interest rates is very unrealistic whereby we are going to approach raise the problem from a different point of view.In effect, the package "Annuity Random Interest Rate", hereinafter referred as AnnuityRIR (Cruz Rambaud et al. 2017b), provides a tool which reflects the different methodologies to approximate the expression of the present and the final value of an annuity when the interest rate is random.In order to do it, we will include the approaches showed in solutions provided by Cruz Rambaud et al. (2015), where it is assumed that all non-central moments of the capitalization factor are known, and Cruz Rambaud et al. (2017a), where it is assumed that some non-central moments of the capitalization factor are unknown.Moreover, in this paper, we will propose some new different treatments of annuities with random interest rate by using R statistical software.Obviously, although unitary annuities are analyzed, the formulas proposed may be applied to calculate the value of other annuities with constant amounts where the interest rate is a random variable.In this paper, the different methodologies to approximate the present value of an annuity have been developed under three main categories: -First, it is assumed that all non-central moments of the capitalization factor are known.
Some previous research focuses on obtaining the moment generating function (De Schepper et al. 1992) and some moment values (Zisheng and Yin 2003).Specifically, in order to calculate the value of these annuities, we propose two different expressions by assuming that the random interest rate is either a discrete or continuous variable.-On the other hand, it is assumed that some non-central moments of negative order are known.In this case, the formula of the mathematical expectation of the ratio of two variables derived by Mood et al. (1974) has been employed to estimate the expected value of a certain annuity.-Finally, if only some positive-order moments are known, the formula by Mood et al. and the cubic discounting, as a more accurate procedure than quadratic discounting, have been employed to calculate the expected present value of an annuity.
This paper is organized as follows.Section 2 introduces the nomenclature and clarifies what of the five aforementioned methodologies have been used to calculate the mathematical expression of the present and final expected value of an n-payment unitary annuity, made at the end/beginning of every year (annuity-immediate and annuity-due, respectively).In Sect.2.1, the expected value of the present value is calculated by considering all non-central moments of negative orders of a discrete (2.1.1)or a continuous random variable (2.1.2).On the other hand, in Sect.2.2, these expected values are calculated when some non-central moments of negative orders are unknown.In Sect.2.3, we are going to consider only moments of positive orders in order to apply the formula by Mood et al. (1974) (2.3.1) or the cubic discounting (2.3.2) as a more accurate procedure than quadratic discounting.Finally, Sect.2.4

Background
In this paper, we will consider the interest rate as a random variable which will be represented as X .Therefore, the capitalization factor, 1 + i, is also a random variable represented as U .Obviously, it is verified that U = 1 + X , so the relationship between means and standard deviations of both variables is: Henceforth, when the mean and standard deviation are mentioned we will refer to the random variable U , unless otherwise specified.
Table 1 exhibits both the present and the future value (PV and FV, respectively) of an n-payment annuity, with payments of 1 unit each made at the end/beginning of every year (annuity-immediate/annuity-due), valued at the rate X = U − 1, approximated by using several methodologies (see Cruz Rambaud et al. 2015;Cruz Rambaud et al. 2017a).
Specifically, the final value of immediate and due annuities can be calculated by the folowing methodologies: Mood et al., quadratic discounting, fitted arctan function, normal distribution and beta distribution.However, the present value of immediate and due annuities only have been calculated by the methodologies due to Mood et al., quadratic discounting and fitted arctan function.In this paper, we will introduce some new methodologies to calculate the present value of an annuity (immediate and due).

Calculation using the discrete random variable
In this case, the present value of an annuity-immediate is the following random variable: Thus, its expected value, which can be computed with the function PV_post_exact(data, years) where μ −r = E(U −r ) is the moment of order −r with respect to the origin of the random variable U ; it can be computed with the function moment(x, order, central, absolute, na.rm).
Observe that, as previously indicated, the non-central moments and so the present value are determined by the information on interest rates at present moment which conditions all subsequent formulation.Therefore, in this case, we have to calculate the negative-order moments of the interest rates distribution by applying the definition to the real data, that is to say, to the distribution of probability obtained from the information on real interest rates.
It is well known that the moment of this discrete random variable adopts has the following expression: where p i is the probability that the random variable takes the value u i .Finally, the expected present value of an annuity-due (PV_pre_exact(data,years)) is : Example 1 Let us consider the annual values of the EURIBOR offered by the main twenty banking institutions of the European Union, on July 31, 2017.By considering a spread of 2%, the interest rate is (see Table 2): In this case, the expected present value of an 10-payment unitary annuity-immediate, valued at the rate X is, by applying Eq. ( 4), whilst the present value of the corresponding annuity-due is, by applying Eq. ( 6), The results can be easily obtained with the examples of the two functions PV_post_exact(data,10) and PV_pre_exact(data,10) which can be found in the manual of the R package "AnnuityRIR ".

Calculation using a continuous random variable
An alternative way to approximate the present value of an annuity is to fit the real interest rates to a continuous distribution function and then to use all the information implicit in the fitted adjusted distribution to calculate the negative-order moments.It is well known that, in the continuous case, the expression of the moment of order −r is: for all values of r , being f (u) the density function of the random variable U .
Example 2 Using the function triangular_parameters_U(data) (Cruz Rambaud et al. 2017b), the data of the interest rates shown in Table 2 can be fitted to a triangular distribution of parameters a = 0.018182, b = 0.018856 and c = 0.018442.Therefore, by applying Eq. ( 7) (function triangular_moments_dis_U(data,order)), the expected present value of an a 10-payment unitary annuity-immediate, valued at the rate X and computed with the function PV_post_triang_dis(data,10) is E(a 10 X ) = 9.053807, whilst the present value of the corresponding annuity-due computed with the function PV_pre_triang_dis(data,10) is Alternatively and like the positive-order moments, we can use another formula to find the negative-order moments of a continuous random variable, according to the generating moments function M U (t) (Meng 2005): Example 3 It can be shown that the interest rates of Table 2 fit a normal distribution of mean μ = 0.01848 and standard deviation σ = 0.00028 (function norm_test_jb(data)).The problem exhibited by the normal distribution is that all its negative-order moments are divergent.Alternatively, we can fit the data to a beta distribution (function beta_parameters(data)) but the calculation of these moments involves solving the convergence of an alternate series.
We can also use the triangular distribution of Example 2, in which the moments of order −1 and −2 are infinite but the rest of negative-order moments (r > 2) can be determined with the function triangular_moments_3_U(data,order).In effect, the moment generating function of the triangular distribution of parameters a (pessimistic), b (optimistic) and c (most likely) is the following: In this case, taking into account that (k > 0): it can be easily shown that (r > 2): In this case, the expected present value of a 10-payment unitary annuity-immediate, valued at the rate X and computed with the function denoted by PV_post_triang_3(data,10) is E(a 10 X ) = 9.053807, whilst the present value of the corresponding annuity-due computed with the function PV_pre_triang_3(data,10) is Both these results can be easily reproduced with the examples which are displayed in the manual of the package "AnnuityRIR" (Cruz Rambaud et al. 2017b).

Considering some non-central moments of negative orders
In this Subsection, we are going to approximate the expression of the present expected value of an n-payment annuity, with payments of 1 unit each made at the end of every year (annuityimmediate), valued at the rate X , by using the formula of the mathematical expectation of the ratio of two random variables X and Y derived by Mood et al. (1974) and introduced as well by Rice (1995): (8) Thus, we can obtain the present value as follows: Thus, taking into account that (Fisz 1963): Analogously, the present expected value of an n-payment annuity, with payments of 1 unit each made at the beginning of every year (annuity-due), valued at the rate X , is: Thus, taking into account that: we can write Example 4 By considering the data of the interest rates shown in Table 2 and by applying Eq. ( 9) with the function PV_post_mood_nm(data,10), the expected present value of a 10-payment unitary annuity-immediate, valued at the rate X , is whilst the present value of the corresponding annuity-due, computed with the function PV_pre_mood_nm(data,10) (Eq.( 10)), is E( ä10 X ) = 9.221975.

Calculation using the Mood et al. formula
In this Subsection, we are going to obtain the expression of the present expected value of an n-payment annuity, with payments of 1 unit each made at the end of every year (annuityimmediate), valued at the rate X , by using now the following expression: which only involves some positive-order moments.As Some algebraic calculation shows that: Analogously, the present expected value of an n-payment annuity, with payments of 1 unit each made at the beginning of every year (annuity-due), valued at the rate X , but can be estimated by using the following expression: Therefore, Some algebraic calculation shows that: Example 5 By considering the data of the interest rates shown in Table 2 and by applying Eq. ( 11) with the function PV_post_mood_pm(data,10), the expected present value of an 10-payment unitary annuity-immediate, valued at the rate X is whilst the present value of the corresponding annuity-due, computed with the function PV_pre_mood_pm(data,10) (Eq.( 12)) is The manual of "AnnuityRIR" (Cruz Rambaud et al. 2017b) allows us to reproduce these results using the examples number 3 of the aforementioned functions.

Calculation using the cubic approximation
In this Subsection, we are going to introduce a new different procedure consisting of the expansion of each power U −r = (1 + X ) −r by using McLaurin's formula: Therefore, Simple algebra shows that: Therefore, Analogously, it can be shown that This approximation exhibits an inconvenience: the expression of the present value is based on the moments of the three first orders and this can distort the result, more specially, when the number of summands is high.Moreover, observe that now the moments involved in ( 14) correspond to the random variable X and not U .
Example 6 By considering again the data of the interest rates shown in Table 2 and by applying Eq. ( 13), the expected present value of a 10-payment unitary annuity-immediate, valued at the rate X is E(a 10 X ) = 9.054472, whilst the present value of the corresponding annuity-due (Eq.( 14)) is These results can be checked with the functions PV_post_cubic(data,10) and PV_pre_cubic(data,10), respectively.
Summarizing, Table 3 shows the results provided by the five used methods.Observe that all them give very similar results in calculating the present value of an 10-payment annuity, with payments of 1 unit each valued at the rate X = U − 1 (with a random interest rate).

An approximation to the variance
Finally, taking into account that the expression of the present value obtained by Eq. ( 1) is an average value, it is necessary to obtain the standard deviation in order to know the possible error in which we have incurred.It is well known that the general expression of the variance of the sum of n random variables X k (independent or not) (15)

Calculation using the discrete random variable
In the particular case in which Y = n k=1 X k , taking into account that var(X k ) = μ 2k − μ 2 k and that cov(X r , X s ) = μ r +s − μ r μ s , where μ k is the non-central moment of order k, one has (n ≤ 2): This formula can be shown by recurrence on n.
Writing the variances and the covariance according to its corresponding non-central moments, Let us suppose that formula (16) holds for n − 1 and let us demonstrate this expression for n.In this case, by using expression (17), one has: Taking into account the recurrence hypothesis and that cov(W + Z , T ) = cov(W , T ) + cov(Z , T ), one has: Observe that μ 2 n and the last summand in the second side complete the square of the sum.Moreover, the penultimate summand completes the second sum.Therefore, this leads to expression (15).
The function variance_drv(data,years) provides the calculation of the variance according to this method.

Calculation using the Mood et al. formula
In this Subsection, we will use the formula of the variance of the ratio of two random variables X and Y derived by Mood et al. (1974):

Conclusion
In the present or final value of an annuity, slight variations of the employed interest rate may involve high variations in the results.Additionally, the estimation of the present and final value of those annuities embedded in uncertain environments, where the employed interest rate is neither constant nor known, entails a high level of difficulty.In this way, this paper takes advantage of the accuracy and the powerful information obtained from the statistical distributions fitted to empirical information on interest rates.
This paper introduces and illustrates some new approaches to estimate the present value of an n-payment unitary annuity with random interest rate.A package of the R statistical software, called "Annuity Random Interest Rate" (hereinafter, AnnuityRIR), has been elaborated as a complementary tool to facilitate the calculation and implementation of these new methodologies.
The approaches proposals presented in this paper have been developed by using different and alternative procedures mainly based on the knowledge of some non-central moments of the random interest rate.More specifically, it should be highlighted the calculation of the expected present value of an annuity starting from the non-central moments of negative order of the random variable (discrete or continuous) describing the interest rate.In effect, negative-order moments have been scarcely investigated and mainly from a theoretical point of view.It is noteworthy to highlight that all involved moments have been calculated by taking the present moment as a benchmark, that is to say, instant 0 has been the reference point for all estimates.This is logical because the present and final values of an annuity are referenced to the origin of such annuity.This paper presents a practical application to calculate the expected present value of an annuity, being, to the extent of our knowledge, one of the few practical applications of negative-order moments.
On the other hand, the expected values have been calculated by using the formula by Mood et al. when some non-central moments of negative order are given.Moreover, this paper introduces the so-called cubic discounting, as a procedure more accurate than quadratic discounting, which have been employed to calculate the present value of an annuity when only positive moments are given.The non-central moments of all (positive and negative) orders have been derived from the real data exhibited by the EURIBOR at the present moment.If possible, the data have been fitted to a well-known statistical distribution, thus taking advantage of the its accuracy of the distribution.This way we will be able to use all the powerful information contained in such distribution in order to calculate its negative-order moments.
Finally, this paper has provided the expressions of the variance of the present value by considering all the presented approaches.As known, this parameter is completely necessary to obtain the (absolute and relative) error of the approximate values obtained in this paper.Clearly, this work is not without limitations; indeed, other distributions may be considered and other non-parametric approaches could be used in future studies.
Funding Open access funding provided by Università degli Studi della Campania Luigi Vanvitelli within the CRUI-CARE Agreement.This paper has been partially supported by the project "La sostenibilidad del Functions Description plot_FVs_pre Plot the final expected values of an n-payment annuity, with payments of 1 unit each made at the beginning of every year (annuity-due), valued at the rate X , using different approaches plot_FV_post_beta_kmom Plot the final expected value of an n-payment annuity, with payments of 1 unit each made at the end of every year (annuity-immediate), valued at the rate X , using the estimated moments of the beta distribution plot_FV_post_norm_kmom Plot the final expected value of an n-payment annuity, with payments of 1 unit each made at the end of every year (annuity-immediate), valued at the rate X , using the estimated moments of the normal distribution plot_FV_pre_beta_kmom Plot the final expected value of an n-payment annuity, with payments of 1 unit each made at the beginning of every year (annuity-due), valued at the rate X , using the estimated moments of the beta distribution plot_FV_pre_norm_kmom Plot the final expected value of an n-payment annuity, with payments of 1 unit each made at the beginning of every year (annuity-due), valued at the rate X , using the estimated moments of the normal distribution plot_PVs_post Plot the present expected values of an n-payment annuity, with payments of 1 unit each made at the end of every year (annuity-immediate), valued at the rate X , using different approaches plot_PVs_pre Plot the present expected values of an n-payment annuity, with payments of 1 unit each made at the beginning of every year (annuity-due), valued at the rate X , using different approaches PV_post_artan Compute present expected value of an n-payment annuity, with payments of 1 unit each, made at the end of every year (annuity-immediate), valued at the rate X , using the tetraparametric function approach PV_post_cubic Compute the present expected value of an n-payment annuity, with payments of 1 unit each made at the end of every year (annuity-due), valued at the rate X , using the cubic discount method PV_post_exact Computes the present value of an annuity-immediate considering only non-central moments of negative orders PV_post_mood_nm Compute the present expected value of an n-payment annuity, with payments of 1 unit each made at the end of every year (annuity-immediate), valued at the rate X , with the method of Mood et al. using some negative moments of the distribution Functions Description triangular_moments_dis Compute the negative moments of the fitted triangular distribution of the random variable X according to the definition (as integral) triangular_moments_dis_U Compute the negative moments of the fitted triangular distribution of the random variable "capitalization factor" U according to the definition (as integral) triangular_parameters Compute the parameters and plot the fitted triangular distribution of the random variable X triangular_parameters_U Return the parameters of the fitted triangular distribution of the random variable "capitalization factor" U variance_drv Compute the variance of the present value of an annuity using "discrete random variable" approach variance_post_mood_nm Compute the variance of the present value of an annuity-immediate using the

Fig. 1
Fig. 1 A summary of presented methodologies.Source: own elaboration

Table 1
Updating the status of the research

Table 2
Data values of EURIBOR on July 31, 2017

Table 3
Summary of the obtained values of the present value Mood et al. approximation and some non-central moments of negative order variance_post_mood_pm Compute the variance of the present value of an annuity-immediate using the Mood et al. approximation and some non-central moments of positive order variance_pre_mood_nm Compute the variance of the present value of an annuity-due using the Mood et al. approximation and some non-central moments of negative order variance_pre_mood_pm Compute the variance of the present value of an annuity-due using the Mood et al. approximation and some non-central moments of positive order