# Pension Funding with lifecontingencies

#### Gian Paolo Clemente, Catholic University, Milan

#### 05 marzo, 2019

## Introduction and some general settings

This report focuses on several funding and actuarial costing methods in occupational pension funds where the benefits are defined. As well-known, a defined benefit (DB) pension plan is a type of pension plan in which an employer/sponsor promises a specified pension payment, lump-sum (or combination thereof) on retirement that is predetermined by a formula based on the employee’s history (earnings, tenure of service and age, etc.) rather than depending directly on individual investment returns. The first consideration is what constitutes the “retirement provision”. We could argue that provisions regard any savings or investments made by an individual for later life. For the sake of brevity, we only focus here on the actuarial involvement in the provision of retirement benefits by employers, groups of employers and professions for their members (i.e the so called “second pillar” of retirement provisions). This provision is usually designed to supplement and in some cases replace (“contracted out” pension schemes in the United Kingdom) the compulsory state pension arrangements. The main benefits are usually at least both a pension payable for life from the attainment of a specified retirement age and a pension payable to the widow/er on the death of the employee/member, according to the death of either the active member or the retired member. These core benefits are found in almost every pension funds, albeit in differing forms. Also additional benefits are usually provided in variable frameworks (e.g. disability for instance).
In the next, we focus on a specific defined benefit pension fund by considering only the retirement benefit for the sake of simplicity. We present some well-known funding methods and we show how to evaluate both premium rate and technical reserve via life contingencies package dedicated functions. To this aim we follow a cohort approach by considering a group of people who share the same characteristic in terms of age of affiliation \(\alpha\) with the pension plan, sex, age of retirement \(\beta\) and initial salary (\(s_\alpha\)). We consider a DB scheme where the pension on retirement is fixed in advance as a proportion of the member’s salary in their last year of service. So the level of the pension at the time of benefit entitlement is based on a proportional-earnings related formula and benefit is expressed as a pension annuity \(b_{x}\), starting at age \(\beta\), increasing annualy at the revaualion rate \(\delta\) and payable until death: \[\begin{equation} b_{x}=(1+\delta)^{x-\beta}\left(\frac{\beta-\alpha}{t}\right)s_{\beta} \end{equation}\] In particular, we are assuming that:

In other words, the pension plan guarantees a replacement ratio equal to: \[\begin{equation} RR_{\beta}=\frac{s_{\beta}}{b_{\beta}}=\left(\frac{\beta-\alpha}{t}\right) \end{equation}\]

Furthermore, we assume that annual earnings are increasing at a yearly flat rate \(i\): \[\begin{equation} s_{x}=s_{x-1}(1+j) \end{equation}\] This evolution of earnings, desibed by the rate \(j\), is usually considered to include two components:

## Contribution rate and technical reserves

Our first aim is to quantify the level of contribution to be set at a level to produce the targeted pension on retirement. In this framework, two general approaches can be pursued. On one hand, the accrued benefits cost method defines the annual normal cost on the basis of portion of pension benefit matured at the valuation date. Accrued Benefits Funding Methods are indeed a major category of funding methods in which the Actuarial Liability for active members is based on pensionable service accrued up to the valuation date or to the end of the Control Period, asappropriate. Differences between the various Accrued Benefits Funding Methods arise from the treatment of decrements in membership and increases in pensionable pay when calculating the Actuarial Liabilities for active members. This affects the value placed not only on the Actuarial Liability but also on the Standard Contribution Rate. Main examples treated in the next are Current Unit Method and Projected Unit Method. \ On the other hand the Projected Benefit Cost Methods, that include several methods (as Individual Entry-Age Method for instance), quantifies a normal cost that represents an amount that will provide for the estimated projected retirement benefits over the service lives of either the individual employees or the employee group. Pension costs computed under this approach tend to be stable or decline year by year, depending on the method selected.

We begin by the Current Unit Method (CUM). The Actuarial Liability for active members is calculated taking into account all types of decrement (in our case only retirement). In calculating the Actuarial Liability as at the valuation date pensionable pay is not projected. It is now possible to define the following mathematical reserve V_{x} at age \(x\) (\(\alpha \leq x <\beta\)): \[\begin{equation} V_{x}={}_{\beta-x}E^{p^{ap},i^{ap}}_{x}\left(\frac{x-\alpha}{t}s_{x-1}\right)\ddot{a}^{p^{pp},u}_{\beta} \end{equation}\] where \(p^{ap}\) and \(p^{pp}\) are the survival probabilities in accumulation and payment periods respectively, \(i^{ap}\) is the technical financial rate in the accumulation period, \(u=\frac{1+i^{pp}}{1+\delta}-1\) is a syntethic index function of the technical financial rate \(i^{pp}\) in the payment period and the indexation of the pension amount. Notice that \(p^{pp}\) and \(j^{pp}\) may depend on premium rate requested by the insurer when an agreement is made between pension plan and insurer for the payments of annuity amounts. \({}_{\beta-x}E^{p^{ap},i^{ap}}_{x}=\frac{{}_{\beta-x}p^{ap}_{x}}{(1+i)^{\beta-x}}\) is the actuarial present value of a life-contingent \(\beta-x\)-year future payment of 1. While \(\ddot{a}^{p^{pp},u}_{\beta}\) is the expected present value for the life-annuity due of infinite-duration. By applying Fouret’s formula, it is easily derive the contribution rate \(\gamma_{x}\) as: \[\begin{equation} \gamma_{x}={}_{\beta-x}E^{p^{ap},i^{ap}}_{x}\frac{1}{t}\left(1+\frac{j}{1+j}(x-\alpha)\right)\ddot{a}^{p^{pp},u}_{\beta} \end{equation}\] In such calculations, allowance is made for increases in the benefits between the relevant date and the assumed date of retirement because of the future increases in salaries. Finally

The Projected Unit Method is based on a similar structure of CUM, but the pensionable pay is projected from the relevant date up to the assumed date of retirement. This method is also known as the Projected Unit Credit Method. We have then the following relations for the mathematical reserve and the contribution rate respectively: \[\begin{equation} V_{x}={}_{\beta-x}E^{p^{ap},i^{ap}}_{x}\left(\frac{x-\alpha}{t}s_{\beta-1}\right)\ddot{a}^{p^{pp},u}_{\beta} \end{equation}\] \[\begin{equation} \gamma_{x}={}_{\beta-x}E^{p^{ap},i^{ap}}_{x}\frac{1}{t}(1+j)^{\beta-x-1}\ddot{a}^{p^{pp},u}_{\beta} \end{equation}\] In this case, the normal cost depends yearly on a portion \(\frac{1}{t}\) of the expected earning at the last year of service before than retirement.

Finally, we consider Individual Entry-Age Method. This method assumes that every employee entered the plan (thus, entry age) at the time of employment and that contributions have been made on this basis from the entry age to the date of valuation. The contributions are the level annual amounts which, if accumulated at the rate of interest used in the actuarial valuation, would result in a fund equal to the present value of the pension benefits at retirement for employees that survive at that time. In this case the contribution rate is equal to the ratio of expected present value of future benefits to expected present value of future earnings: \[\begin{equation} \gamma_{x}=\frac{{}_{\beta-\alpha}E^{p^{ap},i^{ap}}_{\alpha}\left(\frac{\beta-\alpha}{t}s_{\beta-1}\right)\ddot{a}^{p^{pp},u}_{\beta}}{s_{\alpha}{}_{\beta-\alpha}\ddot{a}^{p^{pp},r}_{\alpha}} \end{equation}\] where \(u=\frac{1+i^{pp}}{1+j}-1\) is a syntethic index function of the technical financial rate \(i^{pp}\) in the payment period and the annual salary increase \(j\). In this case, \(\gamma_{x}\) is constant over time.

Mathematical reserve can be derived in a perspective view as: \[\begin{equation} V_{x}={}_{\beta-x}E^{p^{ap},i^{ap}}_{x}\left(\frac{\beta-\alpha}{t}s_{\beta-1}\right)\ddot{a}^{p^{pp},u}_{\beta}-\gamma_{x}s_{x}{}_{\beta-x}\ddot{a}^{p^{pp},r}_{x} \end{equation}\]

## Functions definitions

The CUM contribution rate is defined by following function:

    CUM <- function(acttableAccPeriod, x, beta, i = actuarialtable@interest, j,
        t, k = 1, payment = "advance", acttablePaymPeriod, i2, delta = 0) {

        out <- numeric(1)
        if (missing(acttableAccPeriod))
            stop("Error! Need an actuarial actuarialtable")
        if (missing(acttablePaymPeriod))
            acttablePaymPeriod = acttableAccPeriod
        if (missing(i2))
            i2 = i
        if (missing(x))
            stop("Error! Need age!")
        if (missing(beta))
            stop("Error! Retirement age!")
        if (x > getOmega(acttableAccPeriod)) {
            out = 0
            return(out)
        }
        if (missing(t))
            stop("Error! Need t")
        if (missing(j))
            stop("Error! Need average salary increase rate")
        if (any(x < 0, beta < 0, t < 0))
            stop("Error! Negative parameters")
        out = sapply(seq(x, beta - 1, 1), function(h) Exn(acttableAccPeriod, h,
            beta - h, i = i) * ((1/t) + (h - x)/t * (j/(1 + j))) * axn(acttablePaymPeriod,
            beta, i = (1 + i2)/(1 + delta) - 1, k = 1))
        return(out)
    }

The mathematical reserve is instead equal to:

    CUMmr <- function(acttableAccPeriod, x, beta, i = actuarialtable@interest, j,
        t, k = 1, payment = "advance", acttablePaymPeriod, i2, delta = 0) {

        out <- numeric(1)
        if (missing(acttableAccPeriod))
            stop("Error! Need an actuarial actuarialtable")
        if (missing(acttablePaymPeriod))
            acttablePaymPeriod = acttableAccPeriod
        if (missing(i2))
            i2 = i
        if (missing(x))
            stop("Error! Need age!")
        if (missing(beta))
            stop("Error! Retirement age!")
        if (x > getOmega(acttableAccPeriod)) {
            out = 0
            return(out)
        }
        if (missing(t))
            stop("Error! Need t")
        if (missing(j))
            stop("Error! Need average salary increase rate")
        if (any(x < 0, beta < 0, t < 0))
            stop("Error! Negative parameters")
        out = c(sapply(seq(x, beta, 1), function(h) Exn(acttableAccPeriod, h, beta -
            h, i = i) * ((h - x)/t * (1 + j)^(h - x - 1)) * axn(acttablePaymPeriod,
            beta, i = (1 + i2)/(1 + delta) - 1, k = 1)), sapply(seq(beta + 1, getOmega(acttablePaymPeriod) +
            1, 1), function(h) ((beta - x)/t * (1 + j)^(beta - x - 1)) * (1 + delta)^(h -
            beta) * axn(acttablePaymPeriod, h, i = (1 + i2)/(1 + delta) - 1, k = 1)))
        return(out)
    }

The PUM contribution rate is defined by the following function:

    # Projected Unit Method
    PUM <- function(acttableAccPeriod, x, beta, i = actuarialtable@interest, j,
        t, k = 1, payment = "advance", acttablePaymPeriod, i2, delta = 0) {

        out <- numeric(1)
        if (missing(acttableAccPeriod))
            stop("Error! Need an actuarial actuarialtable")
        if (missing(acttablePaymPeriod))
            acttablePaymPeriod = acttableAccPeriod
        if (missing(i2))
            i2 = i
        if (missing(x))
            stop("Error! Need age!")
        if (missing(beta))
            stop("Error! Retirement age!")
        if (x > getOmega(acttableAccPeriod)) {
            out = 0
            return(out)
        }
        if (missing(t))
            stop("Error! Need t")
        if (missing(j))
            stop("Error! Need average salary increase rate")
        if (any(x < 0, beta < 0, t < 0))
            stop("Error! Negative parameters")
        out = sapply(seq(x, beta - 1, 1), function(h) Exn(acttableAccPeriod, h,
            beta - h, i = i) * 1/t * axn(acttablePaymPeriod, beta, i = (1 + i2)/(1 +
            delta) - 1, k = 1) * (1 + j)^(beta - h - 1))
        return(out)
    }

while the mathematical reserve is

    PUMmr <- function(acttableAccPeriod, x, beta, i = actuarialtable@interest, j,
        t, k = 1, payment = "advance", acttablePaymPeriod, i2, delta = 0) {

        out <- numeric(1)
        if (missing(acttableAccPeriod))
            stop("Error! Need an actuarial actuarialtable")
        if (missing(acttablePaymPeriod))
            acttablePaymPeriod = acttableAccPeriod
        if (missing(i2))
            i2 = i
        if (missing(x))
            stop("Error! Need age!")
        if (missing(beta))
            stop("Error! Retirement age!")
        if (x > getOmega(acttableAccPeriod)) {
            out = 0
            return(out)
        }
        if (missing(t))
            stop("Error! Need t")
        if (missing(j))
            stop("Error! Need average salary increase rate")
        if (any(x < 0, beta < 0, t < 0))
            stop("Error! Negative parameters")
        out = c(sapply(seq(x, beta, 1), function(h) Exn(acttableAccPeriod, h, beta -
            h, i = i) * ((h - x)/t * (1 + j)^(beta - x - 1)) * axn(acttablePaymPeriod,
            beta, i = (1 + i2)/(1 + delta) - 1, k = 1)), sapply(seq(beta + 1, getOmega(acttablePaymPeriod) +
            1, 1), function(h) ((beta - x)/t * (1 + j)^(beta - x - 1)) * (1 + delta)^(h -
            beta) * axn(acttablePaymPeriod, h, i = (1 + i2)/(1 + delta) - 1, k = 1)))
        return(out)
    }

The IEAM is defined as follows:

    # Individual Entry-Age Unit Method Type: 0 constant contribution rate, 1 #
    # Constant Contribution amount (Default is 0)
    IEAM <- function(acttableAccPeriod, x, beta, i = actuarialtable@interest, j,
        t, k = 1, payment = "advance", acttablePaymPeriod, i2, delta = 0, type = 0) {

        out <- numeric(1)
        if (missing(acttableAccPeriod))
            stop("Error! Need an actuarial actuarialtable")
        if (missing(acttablePaymPeriod))
            acttablePaymPeriod = acttableAccPeriod
        if (missing(i2))
            i2 = i
        if (missing(x))
            stop("Error! Need age!")
        if (missing(beta))
            stop("Error! Retirement age!")
        if (x > getOmega(acttableAccPeriod)) {
            out = 0
            return(out)
        }
        if (missing(t))
            stop("Error! Need t")
        if (missing(j))
            stop("Error! Need average salary increase rate")
        if (any(x < 0, beta < 0, t < 0))
            stop("Error! Negative parameters")
        if (type == 0) {
            out = (Exn(acttableAccPeriod, x, beta - x, i = i) * (beta - x)/t * axn(acttablePaymPeriod,
                beta, i = (1 + i2)/(1 + delta) - 1, k = 1) * (1 + j)^(beta - x -
                1))/(axn(acttablePaymPeriod, x, beta - x, i = (1 + i)/(1 + j) -
                1, k = 1, payment = "advance"))
        } else {
            out = ((Exn(acttableAccPeriod, x, beta - x, i = i) * (beta - x)/t *
                axn(acttablePaymPeriod, beta, i = (1 + i2)/(1 + delta) - 1, k = 1) *
                (1 + j)^(beta - x - 1))/(axn(acttablePaymPeriod, x, beta - x, i,
                k = 1, payment = "advance")))/(1 + j)^seq(0, beta - x - 1, 1)
        }
        return(out)
    }

while the reserve is:

    IEAMmr <- function(acttableAccPeriod, x, beta, i = actuarialtable@interest,
        j, t, k = 1, payment = "advance", acttablePaymPeriod, i2, delta = 0, type = 0) {

        out <- numeric(1)
        if (missing(acttableAccPeriod))
            stop("Error! Need an actuarial actuarialtable")
        if (missing(acttablePaymPeriod))
            acttablePaymPeriod = acttableAccPeriod
        if (missing(i2))
            i2 = i
        if (missing(x))
            stop("Error! Need age!")
        if (missing(beta))
            stop("Error! Retirement age!")
        if (x > getOmega(acttableAccPeriod)) {
            out = 0
            return(out)
        }
        if (missing(t))
            stop("Error! Need t")
        if (missing(j))
            stop("Error! Need average salary increase rate")
        if (any(x < 0, beta < 0, t < 0))
            stop("Error! Negative parameters")
        al = IEAM(acttableAccPeriod, x, beta, i, j, t, k = 1, payment, acttablePaymPeriod,
            i2, delta, type)
        if (type == 0) {
            out = c(sapply(seq(x, beta, 1), function(h) Exn(acttableAccPeriod, h,
                beta - h, i = i) * ((beta - x)/t * (1 + j)^(beta - x - 1)) * axn(acttablePaymPeriod,
                beta, i = (1 + i2)/(1 + delta) - 1, k = 1) - al * (1 + j)^(h - x) *
                axn(acttableAccPeriod, h, beta - h, i = (1 + i)/(1 + j) - 1, k = 1)),
                sapply(seq(beta + 1, getOmega(acttablePaymPeriod) + 1, 1), function(h) ((beta -
                    x)/t * (1 + j)^(beta - x - 1)) * (1 + delta)^(h - beta) * axn(acttablePaymPeriod,
                    h, i = (1 + i2)/(1 + delta) - 1, k = 1)))
        } else {
            out = c(sapply(seq(x, beta - 1, 1), function(h) Exn(acttableAccPeriod,
                h, beta - h, i = i) * ((beta - x)/t * (1 + j)^(beta - x - 1)) *
                axn(acttablePaymPeriod, beta, i = (1 + i2)/(1 + delta) - 1, k = 1) -
                al[h - x + 1] * (1 + j)^(h - x) * axn(acttableAccPeriod, h, beta -
                    h, i, k = 1)), sapply(seq(beta, getOmega(acttablePaymPeriod) +
                1, 1), function(h) ((beta - x)/t * (1 + j)^(beta - x - 1)) * (1 +
                delta)^(h - beta) * axn(acttablePaymPeriod, h, i = (1 + i2)/(1 +
                delta) - 1, k = 1)))
        }
        return(out)
    }

## An applied example

We consider a pension fund based on a single cohort of age \(x=25\). Furthermore, we assume that the DB pension fund has based the quantification of the contribution rate on the following assumptions:

    #Current Unit Method
    beta=65 # Beta Retirement age
    x=25 #x Age of the insured.
    i=0.08 # Interest Rate
    t=60 #1/t is the % of the salary, recognized as retirement pension, for each year of service
    j=0.06 #  average salary increases (for both growth in wages and promotional salary for seniority)
    delta=0.03 #Increase of retirement pension

Therefore, the following Figure shows the pattern of the contribution rates according to the three methods. As expected, Current Unit Method, being based on the current salary, shows a very increasing tendency over time. Differences with respect to PUM depends on the value of \(j\). It is indeed easy to prove that both methods lead to the same rates when \(j=0\).

    CUM(lt, x, beta, i, j, t, k, delta = 0.03)

    ##  [1] 0.005754777 0.006573606 0.007487168 0.008505653 0.009639970
    ##  [6] 0.010902408 0.012306188 0.013866649 0.015599806 0.017523147
    ## [11] 0.019657605 0.022024425 0.024649058 0.027559785 0.030787341
    ## [16] 0.034364816 0.038330247 0.042724286 0.047593736 0.052990323
    ## [21] 0.058977977 0.065625324 0.073010784 0.081220029 0.090346612
    ## [26] 0.100492845 0.111794330 0.124400086 0.138480677 0.154222894
    ## [31] 0.171849693 0.191613090 0.213819124 0.238777442 0.266849269
    ## [36] 0.298627177 0.334515275 0.375154465 0.421219575 0.473518531

    PUM(lt, x, beta, i, j, t, k, delta = 0.03)

    ##  [1] 0.05584152 0.05695271 0.05808426 0.05923836 0.06041482 0.06161534
    ##  [7] 0.06283977 0.06409193 0.06537115 0.06667668 0.06801402 0.06938191
    ## [13] 0.07078554 0.07222979 0.07371758 0.07524967 0.07683003 0.07845965
    ## [19] 0.08014290 0.08188362 0.08369485 0.08558446 0.08756199 0.08963405
    ## [25] 0.09180487 0.09407693 0.09647185 0.09900614 0.10169720 0.10455753
    ## [31] 0.10760721 0.11086517 0.11436081 0.11810296 0.12210671 0.12646536
    ## [37] 0.13115460 0.13622408 0.14170129 0.14762637

    IEAM(lt, x, beta, i, j, t, k, delta = 0.03, type = 0)

    ## [1] 0.08220932

    plot(seq(x, beta - 1, 1), CUM(lt, x, beta, i, j, t, k, delta = 0.03), xlab = "age",
        ylab = "Contribution rate")
    lines(seq(x, beta - 1, 1), PUM(lt, x, beta, i, j, t, k, delta = 0.03), type = "p",
        col = "red")
    lines(seq(x, beta - 1, 1), rep(IEAM(lt, x, beta, i, j, t, k, delta = 0.03),
        beta - x), type = "p", col = "blue")
    lines(seq(x, beta - 1, 1), IEAM(lt, x, beta, i, j, t, k, delta = 0.03, type = 1),
        type = "p", col = "green")
    legend("topleft", c("CUM", "PUM", "IEAM (constant rate)", "IEAM (constant premium)"),
        col = c("black", "red", "blue", "green"), pch = c(1, 1, 1, 1), cex = 0.6)

According to the mathematical reserve, we observe in the following figure the different behaviour in the accumulation period.

    CUMmr(lt, x, beta, i, j, t, k, delta = 0.03)

    ##  [1]  0.000000000  0.006221448  0.014258614  0.024508958  0.037446891
    ##  [6]  0.053639272  0.073760123  0.098616348  0.129162182  0.166528061
    ## [11]  0.212070940  0.267382722  0.334373579  0.415314980  0.512895315
    ## [16]  0.630285313  0.771265800  0.940286298  1.142652316  1.384646537
    ## [21]  1.673897013  2.019416688  2.431989936  2.924395137  3.511752274
    ## [26]  4.211936519  5.047128122  6.043771192  7.233706316  8.654837296
    ## [31] 10.353327248 12.384709961 14.817261640 17.730754052 21.221821606
    ## [36] 25.422343025 30.470089659 36.547276205 43.870093002 52.704951088
    ## [41] 63.404228120 63.168110190 62.847586247 62.445473398 61.946937115
    ## [46] 61.363427358 60.732489682 60.048636888 59.331630821 58.576086731
    ## [51] 57.777042337 56.931141293 56.093837346 55.223520138 54.310320836
    ## [56] 53.357611214 52.375814175 51.403062922 50.413605596 49.413915090
    ## [61] 48.406800653 47.396367590 46.389297076 45.382944194 44.383483293
    ## [66] 43.399243059 42.433645052 41.478016617 40.541725862 39.629586530
    ## [71] 38.734972023 37.893648177 37.079193410 36.256504380 35.359572414
    ## [76] 34.607989694 33.466402838 32.512131687 28.513480477 20.487513359
    ## [81]  0.000000000

    PUMmr(lt, x, beta, i, j, t, k, delta = 0.03)

    ##  [1]  0.00000000  0.06036987  0.13052695  0.21166149  0.30508925
    ##  [6]  0.41227610  0.53483652  0.67459400  0.83353351  1.01383960
    ## [11]  1.21802757  1.44878331  1.70921298  2.00279251  2.33335721
    ## [16]  2.70510308  3.12280482  3.59165827  4.11759038  4.70719397
    ## [21]  5.36841449  6.10994613  6.94172409  7.87473164  8.92108580
    ## [26] 10.09415097 11.41107194 12.89093037 14.55563833 16.42946516
    ## [31] 18.54123225 20.92370689 23.61646393 26.66049831 30.10355955
    ## [36] 34.02082968 38.46778616 43.52839072 49.29243650 55.86724815
    ## [41] 63.40422812 63.16811019 62.84758625 62.44547340 61.94693711
    ## [46] 61.36342736 60.73248968 60.04863689 59.33163082 58.57608673
    ## [51] 57.77704234 56.93114129 56.09383735 55.22352014 54.31032084
    ## [56] 53.35761121 52.37581418 51.40306292 50.41360560 49.41391509
    ## [61] 48.40680065 47.39636759 46.38929708 45.38294419 44.38348329
    ## [66] 43.39924306 42.43364505 41.47801662 40.54172586 39.62958653
    ## [71] 38.73497202 37.89364818 37.07919341 36.25650438 35.35957241
    ## [76] 34.60798969 33.46640284 32.51213169 28.51348048 20.48751336
    ## [81]  0.00000000

    IEAMmr(lt, x, beta, i, j, t, k, delta = 0.03, type = 0)

    ##  [1]  0.00000000  0.08887591  0.19028587  0.30556880  0.43618413
    ##  [6]  0.58374352  0.74999737  0.93691407  1.14659533  1.38132825
    ## [11]  1.64375310  1.93661706  2.26311151  2.62676757  3.03142210
    ## [16]  3.48120425  3.98080069  4.53521697  5.15015243  5.83187085
    ## [21]  6.58800961  7.42674015  8.35733396  9.38985722 10.53522457
    ## [26] 11.80528100 13.21567882 14.78347904 16.52807339 18.47048884
    ## [31] 20.63572011 23.05177588 25.75258410 28.77174653 32.14826059
    ## [36] 35.94797752 40.21237306 45.01009766 50.41185014 56.50202292
    ## [41] 63.40422812 63.16811019 62.84758625 62.44547340 61.94693711
    ## [46] 61.36342736 60.73248968 60.04863689 59.33163082 58.57608673
    ## [51] 57.77704234 56.93114129 56.09383735 55.22352014 54.31032084
    ## [56] 53.35761121 52.37581418 51.40306292 50.41360560 49.41391509
    ## [61] 48.40680065 47.39636759 46.38929708 45.38294419 44.38348329
    ## [66] 43.39924306 42.43364505 41.47801662 40.54172586 39.62958653
    ## [71] 38.73497202 37.89364818 37.07919341 36.25650438 35.35957241
    ## [76] 34.60798969 33.46640284 32.51213169 28.51348048 20.48751336
    ## [81]  0.00000000

    IEAMmr(lt, x, beta, i, j, t, k, delta = 0.03, type = 1)

    ##  [1]  0.0000000  0.1913544  0.3982143  0.6218430  0.8635912  1.1249462
    ##  [7]  1.4074889  1.7130265  2.0434158  2.4006461  2.7871081  3.2051458
    ## [13]  3.6576036  4.1476131  4.6785141  5.2538136  5.8775626  6.5539643
    ## [19]  7.2878914  8.0846533  8.9510620  9.8941995 10.9221420 12.0434996
    ## [25] 13.2674524 14.6038196 16.0665212 17.6704499 19.4325108 21.3707116
    ## [31] 23.5067767 25.8648860 28.4747820 31.3646959 34.5676347 38.1439516
    ## [37] 42.1263791 46.5748283 51.5494943 57.1227408 63.4042281 63.1681102
    ## [43] 62.8475862 62.4454734 61.9469371 61.3634274 60.7324897 60.0486369
    ## [49] 59.3316308 58.5760867 57.7770423 56.9311413 56.0938373 55.2235201
    ## [55] 54.3103208 53.3576112 52.3758142 51.4030629 50.4136056 49.4139151
    ## [61] 48.4068007 47.3963676 46.3892971 45.3829442 44.3834833 43.3992431
    ## [67] 42.4336451 41.4780166 40.5417259 39.6295865 38.7349720 37.8936482
    ## [73] 37.0791934 36.2565044 35.3595724 34.6079897 33.4664028 32.5121317
    ## [79] 28.5134805 20.4875134  0.0000000

    plot(seq(x, getOmega(lt) + 1, 1), CUMmr(lt, x, beta, i, j, t, k, delta = 0.03),
        xlab = "age", ylab = "", main = "Mathematical reserve")
    lines(seq(x, getOmega(lt) + 1, 1), PUMmr(lt, x, beta, i, j, t, k, delta = 0.03),
        type = "p", col = "red")
    lines(seq(x, getOmega(lt) + 1, 1), IEAMmr(lt, x, beta, i, j, t, k, delta = 0.03,
        type = 0), type = "p", col = "blue")
    lines(seq(x, getOmega(lt) + 1, 1), IEAMmr(lt, x, beta, i, j, t, k, delta = 0.03,
        type = 1), type = "p", col = "green")
    legend("topleft", c("CUM", "PUM", "IEAM (constant rate)", "IEAM (constant premium)"),
        col = c("black", "red", "blue", "green"), pch = c(1, 1, 1, 1), cex = 0.6)