diff --git a/NAMESPACE b/NAMESPACE
index 497e9a369bf72f481f36c72d14bb39173467f929..2a20d65225811de469c0c720bdab9a40687b41e4 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -51,6 +51,7 @@ export(PP.rate.terminalBonus)
export(PP.rate.terminalBonusFund)
export(PP.rate.totalInterest)
export(PP.rate.totalInterest2)
+export(PVfactory)
export(PaymentTimeEnum)
export(ProfitComponentsEnum)
export(ProfitParticipation)
diff --git a/R/HelperFunctions.R b/R/HelperFunctions.R
index 4f68bb23194392768dac4759acdac36032ff3538..ef96a531ed5a3ae0e0348f7fd1752bf5e7dc23b4 100644
--- a/R/HelperFunctions.R
+++ b/R/HelperFunctions.R
@@ -278,185 +278,348 @@ mergeValues3D = function(starting, ending, t) {
abind::abind(starting[1:t,,,], ending[-1:-t,,,], along = 1)
}
}
-# Caution: px is not neccessarily 1-qx, because we might also have dread diseases so that px=1-qx-ix! However, the ix is not used for the survival present value
-calculatePVSurvival = function(px = 1 - qx, qx = 1 - px, advance = NULL, arrears = NULL, ..., m = 1, mCorrection = list(alpha = 1, beta = 0), v = 1, start = 0) {
- init = advance[1]*0;
- if (!any(advance != 0 ) && !any(arrears != 0)) {
- return(init);
- }
- # assuming advance and arrears have the same dimensions...
- l = max(length(qx), dim(advance)[[1]], dim(arrears)[[1]]);
- p = pad0(px, l, value=0);
- if (missing(advance)) {
- advance = arrears * 0
- }
- if (missing(arrears)) {
- arrears = advance * 0
- }
- advance = pad0(advance, l, value = init);
- arrears = pad0(arrears, l, value = init);
-
- # TODO: Make this work for matrices (i.e. currently advance and arrears are assumed to be one-dimensional vectors)
- # TODO: Replace loop by better way (using Reduce?)
- res = rep(0, l + 1);
- advcoeff = mCorrection$alpha - mCorrection$beta*(1-p*v);
- arrcoeff = mCorrection$alpha - (mCorrection$beta + 1/m)*(1-p*v);
- for (i in l:(start + 1)) {
- # coefficients for the payments (including corrections for payments during the year (using the alpha(m) and beta(m)):
- # The actual recursion:
- res[i] = advance[i]*advcoeff[i] + arrears[i]*arrcoeff[i] + v*p[i]*res[i+1];
- }
- res[1:l]
-}
-# Caution: px is not neccessarily 1-qx, because we might also have dread diseases so that px=1-qx-ix! However, the ix is not used for the survival present value
-calculatePVSurvival2D = function(px = 1 - qx, qx = 1 - px, advance = NULL, arrears = NULL, ..., m = 1, mCorrection = list(alpha = 1, beta = 0), v = 1, start = 0) {
-# browser()
- if (!any(advance != 0 ) && !any(arrears != 0)) {
- return(advance * 0);
- }
-
- l = max(length(qx), dim(advance)[[1]], dim(arrears)[[1]]);
- p = pad0(px, l, value=0);
- if (missing(advance)) {
- advance = arrears * 0
- }
- if (missing(arrears)) {
- arrears = advance * 0
- }
- # assuming advance and arrears have the same dimensions...
- advance = padArray(advance, len = l);
- arrears = padArray(arrears, len = l);
-
-
- # TODO: Make this work for arbitrary dimensions
- # TODO: Replace loop by better way (using Reduce?)
- res = advance * 0;
- advcoeff = mCorrection$alpha - mCorrection$beta*(1-p*v);
- arrcoeff = mCorrection$alpha - (mCorrection$beta + 1/m)*(1-p*v);
- for (i in (l-1):start) {
- # coefficients for the payments (including corrections for payments during the year (using the alpha(m) and beta(m)):
- # The actual recursion:
- res[i,] = advance[i,]*advcoeff[i] + arrears[i,]*arrcoeff[i] + v*p[i]*res[i+1,];
- }
- res[1:l,]
-}
-
-
-calculatePVGuaranteed = function(advance, arrears = c(0), ..., m = 1, mCorrection = list(alpha = 1, beta = 0), v = 1, start = 0) {
- # assuming advance and arrears have the same dimensions...
- init = advance[1]*0;
- if (!any(advance != 0 ) && !any(arrears != 0)) {
- return(init);
- }
- l = max(length(advance), length(arrears));
- advance = pad0(advance, l, value = init);
- arrears = pad0(arrears, l, value = init);
-
- # TODO: Make this work for matrices (i.e. currently advance and arrears are assumed to be one-dimensional vectors)
- # TODO: Replace loop by better way (using Reduce?)
- res = rep(0, l + 1);
- advcoeff = mCorrection$alpha - mCorrection$beta * (1 - v);
- arrcoeff = mCorrection$alpha - (mCorrection$beta + 1 / m) * (1 - v);
- for (i in l:(start + 1)) {
- # coefficients for the payments (including corrections for payments during the year (using the alpha(m) and beta(m)):
- # The actual recursion:
- res[i] = advance[i]*advcoeff + arrears[i]*arrcoeff + v*res[i + 1];
- }
- res[1:l]
-}
-
-calculatePVDeath = function(px, qx, benefits, ..., v = 1, start = 0) {
- init = benefits[1]*0; # Preserve the possible array structure of the benefits -> vectorized calculations possible!
- if (!any(benefits != 0 )) {
- return(init);
- }
- l = max(length(qx), length(benefits));
- q = pad0(qx, l, value = 1);
- p = pad0(px, l, value = 0);
- benefits = pad0(benefits, l, value = init);
-
- # TODO: Make this work for matrices (i.e. currently benefits are assumed to be one-dimensional vectors)
- # TODO: Replace loop by better way (using Reduce?)
- res = rep(init, l + 1);
- for (i in l:(start + 1)) {
- # Caution: p_x is not neccessarily 1-q_x, because we might also have dread diseases, so that px=1-qx-ix!
- res[i] = v * q[i] * benefits[i] + v * p[i] * res[i + 1];
- }
- res[1:l]
-}
-
-calculatePVDisease = function(px = 1 - qx - ix, qx = 1 - ix - px, ix = 1 - px - qx, benefits, ..., v = 1, start = 0) {
- init = benefits[1] * 0;
- if (!any(benefits != 0 )) {
- return(init);
- }
- l = min(length(ix), length(qx), length(benefits));
- qx = pad0(qx, l, value = 1);
- ix = pad0(ix, l, value = 0);
- px = pad0(px, l, value = 0);
- benefits = pad0(benefits, l, value = init);
-
- # TODO: Make this work for matrices (i.e. currently benefits are assumed to be one-dimensional vectors)
- # TODO: Replace loop by better way (using Reduce?)
- res = rep(init, l + 1);
- for (i in l:(start + 1)) {
- res[i] = v * ix[i] * benefits[i] + v * px[i] * res[i + 1];
- }
- res[1:l]
-}
-
-# TODO: So far, we are assuming, the costs array has sufficient time steps and does not need to be padded!
-calculatePVCosts = function(px = 1 - qx, qx = 1 - px, costs, ..., v = 1, start = 0) {
- l = max(length(qx), dim(costs)[1]);
- p = pad0(px, l, value = 0);
- costs = costs[1:l,,,];
-
- # NOTICE: Costs that apply to benefits are a special case and cannot be handled
- # here. Instead, their PV is calculated like the death/survival/invalidity benefits
- # in the tariff's PV calculation function. The values calculated here automatically
- # WILL be ignored
-
- # Take the array structure from the cash flow array and initialize it with 0
- res = costs*0;
-
- # only calculate the loop if we have any relevant cost parameter != 0...
-
- # survival cash flows (until death)
- if (any(costs[,,,"survival"] != 0)) {
- prev = res[1,,,"survival"]*0;
- # Backward recursion starting from the last time:
- # Survival Cash flows
- for (i in l:(start + 1)) {
- # cat("values at iteration ", i, ": ", v, q[i], costs[i,,], prev);
- res[i,,,"survival"] = costs[i,,,"survival"] + v*p[i]*prev;
- prev = res[i,,,"survival"];
- }
- }
- # guaranteed cash flows (even after death)
- if (any(costs[,,,"guaranteed"] != 0)) {
- prev = res[1,,,"guaranteed"]*0;
- for (i in l:(start + 1)) {
- res[i,,,"guaranteed"] = costs[i,,,"guaranteed"] + v*prev;
- prev = res[i,,,"guaranteed"];
- }
- }
- # Cash flows only after death
- # This case is more complicated, as we have two possible states of
- # payments (present value conditional on active, but payments only when
- # dead => need to write the Thiele difference equations as a pair of
- # recursive equations rather than a single recursive formula...)
- if (any(costs[,,,"after.death"] != 0)) {
- prev = res[1,,,"after.death"]*0;
- prev.dead = res[1,,,1]*0
- for (i in l:(start + 1)) {
- res[i,,,"after.death"] = p[i] * v * prev + (1 - p[i]) * v * prev.dead
- prev = res[i,,,"after.death"]
- prev.dead = costs[i,,,"after.death"] + v * prev.dead
+# PVfactory (R6Class for present values with arbitrary dimensions) ####
+#' @export
+PVfactory = R6Class(
+ "PVfactory",
+
+ ######################### PUBLIC METHODS ################################# #
+ public = list(
+ initialize = function(qx, m = 1, mCorrection = list(alpha = 1, beta = 0), v = 1) {
+ private$qx = qx;
+ private$m = m;
+ private$mCorrection = mCorrection;
+ private$v = v;
+ },
+ guaranteed = function(advance = NULL, arrears = NULL, start = 0, ..., m = private$m, mCorrection = private$mCorrection, v = private$v) {
+ # General Note: Since the CF vectors can have an arbitrary number of
+ # dimensions, we cannot directly access them via advance[1,..]. Rather,
+ # we have to construct the `[` function manually as a quoted expression,
+ # inserting the required number of dimensions and then evaluating that
+ # expression. This makes this function a little harder to read, but the
+ # performance should not take a hit and it is implemented in a very
+ # general way.
+
+ cfs = list(advance, arrears)
+ # https://stackoverflow.com/a/16896422/920231
+ deflt = cfs[!unlist(lapply(cfs, is.null))][[1]] * 0
+
+ if (missing(advance) || is.null(advance)) advance = deflt
+ if (missing(arrears) || is.null(arrears)) arrears = deflt
+
+ l = max(unlist(lapply(cfs, function(cf) if(!is.null(dim(cf))) dim(cf)[[1]] else length(cf))))
+
+ # TODO: Make sure all CF tensors have the same number of dimensions
+ dms = if (is.null(dim(advance))) length(advance + 1) else dim(advance);
+
+ # Resulting PV tensor has one timestep more than the CF tensors!
+ dms[1] = dms[1] + 1
+ dmNr = if (is.null(dim(advance))) 1 else length(dim(advance))
+
+ # To be able to access the CF tensors in arbitrary dimensions, we
+ # construct the [..] operator manually by quoting it and then inserting
+ # arguments matching the number of dimensions
+ Qadv = Quote(advance[] )[c(1,2,rep(3, dmNr))];
+ Qarr = Quote(arrears[] )[c(1,2,rep(3, dmNr))];
+ Qres = Quote(res[])[c(1,2,rep(3, dmNr))]; # Access the correct number of dimensions
+ QresAss = Quote(res <- tmp)
+
+ VL = function(quoted, time) {
+ eval(quoted %>% `[<-`(3, time))
+ }
+
+ init = VL(Qadv, 1) * 0;
+
+ # assuming advance and arrears have the same dimensions...
+ # TODO: Pad to length l
+ # advance = pad0(advance, l, value = init);
+ # arrears = pad0(arrears, l, value = init);
+
+ # TODO: Replace loop by better way (using Reduce?)
+ res = array(0, dim = dms)
+
+ # Starting value for the recursion:
+ tmp = init
+ QresAss[[2]] = Qres %>% `[<-`(3, dms[[1]])
+ eval(QresAss)
+
+ advcoeff = mCorrection$alpha - mCorrection$beta * (1 - v);
+ arrcoeff = mCorrection$alpha - (mCorrection$beta + 1/m) * ( 1 - v);
+ for (i in l:(start + 1)) {
+ # coefficients for the payments (including corrections for payments during the year (using the alpha(m) and beta(m)):
+ # The actual recursion:
+ tmp = VL(Qadv, i) * advcoeff + VL(Qarr, i) * arrcoeff + v * VL(Qres, i+1);
+ # Assign tmp to the slice res[i, ....]
+ QresAss[[2]] = Qres %>% `[<-`(3, i)
+ eval(QresAss)
+ }
+ VL(Qres, list(1:l))
+ },
+
+ survival = function(advance = NULL, arrears = NULL, start = 0, ..., m = private$m, mCorrection = private$mCorrection, v = private$v) {
+ # General Note: Since the CF vectors can have an arbitrary number of
+ # dimensions, we cannot directly access them via advance[1,..]. Rather,
+ # we have to construct the `[` function manually as a quoted expression,
+ # inserting the required number of dimensions and then evaluating that
+ # expression. This makes this function a little harder to read, but the
+ # performance should not take a hit and it is implemented in a very
+ # general way.
+
+ cfs = list(advance, arrears)
+ # https://stackoverflow.com/a/16896422/920231
+ deflt = cfs[!unlist(lapply(cfs, is.null))][[1]] * 0
+
+ if (missing(advance) || is.null(advance)) advance = deflt
+ if (missing(arrears) || is.null(arrears)) arrears = deflt
+
+ l = max(unlist(lapply(cfs, function(cf) if(!is.null(dim(cf))) dim(cf)[[1]] else length(cf))))
+
+ # TODO: Make sure all CF tensors have the same number of dimensions
+ dms = if (is.null(dim(advance))) length(advance) else dim(advance);
+
+ # Resulting PV tensor has one timestep more than the CF tensors!
+ dms[1] = dms[1] + 1
+ dmNr = if (is.null(dim(advance))) 1 else length(dim(advance))
+
+ # To be able to access the CF tensors in arbitrary dimensions, we
+ # construct the [..] operator manually by quoting it and then inserting
+ # arguments matching the number of dimensions
+ Qadv = Quote(advance[] )[c(1,2,rep(3, dmNr))];
+ Qarr = Quote(arrears[] )[c(1,2,rep(3, dmNr))];
+ Qres = Quote(res[])[c(1,2,rep(3, dmNr))]; # Access the correct number of dimensions
+ QresAss = Quote(res <- tmp)
+
+ VL = function(quoted, time) {
+ eval(quoted %>% `[<-`(3, time))
+ }
+
+ init = VL(Qadv, 1) * 0;
+
+ # assuming advance and arrears have the same dimensions...
+ p = pad0(private$qx$px, l, value=0);
+ # TODO: Pad to length l
+ # advance = pad0(advance, l, value = init);
+ # arrears = pad0(arrears, l, value = init);
+
+ # TODO: Replace loop by better way (using Reduce?)
+ res = array(0, dim = dms)
+
+ # Starting value for the recursion:
+ tmp = init
+ QresAss[[2]] = Qres %>% `[<-`(3, dms[[1]])
+ eval(QresAss)
+
+ advcoeff = mCorrection$alpha - mCorrection$beta * (1 - p * v);
+ arrcoeff = mCorrection$alpha - (mCorrection$beta + 1/m) * (1 - p * v);
+ for (i in l:(start + 1)) {
+ # coefficients for the payments (including corrections for payments during the year (using the alpha(m) and beta(m)):
+ # The actual recursion:
+ tmp = VL(Qadv, i) * advcoeff[i] + VL(Qarr, i) * arrcoeff[i] + v * p[i] * VL(Qres, i+1);
+ # Assign tmp to the slice res[i, ....]
+ QresAss[[2]] = Qres %>% `[<-`(3, i)
+ eval(QresAss)
+ }
+ VL(Qres, list(1:l))
+ },
+
+ death = function(benefits, start = 0, ..., v = private$v) {
+ # General Note: Since the CF vectors can have an arbitrary number of
+ # dimensions, we cannot directly access them via advance[1,..]. Rather,
+ # we have to construct the `[` function manually as a quoted expression,
+ # inserting the required number of dimensions and then evaluating that
+ # expression. This makes this function a little harder to read, but the
+ # performance should not take a hit and it is implemented in a very
+ # general way.
+
+ cfs = list(benefits)
+ if (missing(benefits) || is.null(benefits)) return(0);
+
+ l = max(unlist(lapply(cfs, function(cf) if(!is.null(dim(cf))) dim(cf)[[1]] else length(cf))))
+
+ # TODO: Make sure all CF tensors have the same number of dimensions
+ dms = if (is.null(dim(benefits))) length(benefits) else dim(benefits);
+
+ # Resulting PV tensor has one timestep more than the CF tensors!
+ dms[1] = dms[1] + 1
+ dmNr = if (is.null(dim(benefits))) 1 else length(dim(benefits))
+
+ # To be able to access the CF tensors in arbitrary dimensions, we
+ # construct the [..] operator manually by quoting it and then inserting
+ # arguments matching the number of dimensions
+ Qben = Quote(benefits[] )[c(1,2,rep(3, dmNr))];
+ Qres = Quote(res[])[c(1,2,rep(3, dmNr))]; # Access the correct number of dimensions
+ QresAss = Quote(res <- tmp)
+
+ VL = function(quoted, time) {
+ eval(quoted %>% `[<-`(3, time))
+ }
+
+ init = VL(Qben, 1) * 0;
+
+ # assuming advance and arrears have the same dimensions...
+ p = pad0(private$qx$px, l, value = 0);
+ q = pad0(private$qx$qx, l, value = 1);
+ # TODO: Pad to length l
+ # benefits = pad0(benefits, l, value = init);
+
+ # TODO: Replace loop by better way (using Reduce?)
+ res = array(0, dim = dms)
+
+ # Starting value for the recursion:
+ tmp = init
+ QresAss[[2]] = Qres %>% `[<-`(3, dms[[1]])
+ eval(QresAss)
+
+ for (i in l:(start + 1)) {
+ # coefficients for the payments (including corrections for payments during the year (using the alpha(m) and beta(m)):
+ # The actual recursion:
+ tmp = v * q[i] * VL(Qben, i) + v * p[i] * VL(Qres, i+1);
+ # Assign tmp to the slice res[i, ....]
+ QresAss[[2]] = Qres %>% `[<-`(3, i)
+ eval(QresAss)
+ }
+ VL(Qres, list(1:l))
+ },
+ disease = function(benefits, start = 0, ..., v = private$v) {
+ # General Note: Since the CF vectors can have an arbitrary number of
+ # dimensions, we cannot directly access them via advance[1,..]. Rather,
+ # we have to construct the `[` function manually as a quoted expression,
+ # inserting the required number of dimensions and then evaluating that
+ # expression. This makes this function a little harder to read, but the
+ # performance should not take a hit and it is implemented in a very
+ # general way.
+
+ cfs = list(benefits)
+ if (missing(benefits) || is.null(benefits)) return(0);
+
+ l = max(unlist(lapply(cfs, function(cf) if(!is.null(dim(cf))) dim(cf)[[1]] else length(cf))))
+
+ # TODO: Make sure all CF tensors have the same number of dimensions
+ dms = if (is.null(dim(benefits))) length(benefits) else dim(benefits);
+
+ # Resulting PV tensor has one timestep more than the CF tensors!
+ dms[1] = dms[1] + 1
+ dmNr = if (is.null(dim(benefits))) 1 else length(dim(benefits))
+
+ # To be able to access the CF tensors in arbitrary dimensions, we
+ # construct the [..] operator manually by quoting it and then inserting
+ # arguments matching the number of dimensions
+ Qben = Quote(benefits[] )[c(1,2,rep(3, dmNr))];
+ Qres = Quote(res[])[c(1,2,rep(3, dmNr))]; # Access the correct number of dimensions
+ QresAss = Quote(res <- tmp)
+
+ VL = function(quoted, time) {
+ eval(quoted %>% `[<-`(3, time))
+ }
+
+ init = VL(Qben, 1) * 0;
+
+ # assuming advance and arrears have the same dimensions...
+ p = pad0(private$qx$px, l, value = 0);
+ ix = pad0(private$qx$ix, l, value = 0);
+ # TODO: Pad to length l
+ # benefits = pad0(benefits, l, value = init);
+
+ # TODO: Replace loop by better way (using Reduce?)
+ res = array(0, dim = dms)
+
+ # Starting value for the recursion:
+ tmp = init
+ QresAss[[2]] = Qres %>% `[<-`(3, dms[[1]])
+ eval(QresAss)
+
+ for (i in l:(start + 1)) {
+ # coefficients for the payments (including corrections for payments during the year (using the alpha(m) and beta(m)):
+ # The actual recursion:
+ tmp = v * ix[i] * VL(Qben, i) + v * p[i] * VL(Qres, i+1);
+ # Assign tmp to the slice res[i, ....]
+ QresAss[[2]] = Qres %>% `[<-`(3, i)
+ eval(QresAss)
+ }
+ VL(Qres, list(1:l))
+ },
+ # Cash flows only after death
+ # This case is more complicated, as we have two possible states of
+ # payments (present value conditional on active, but payments only when
+ # dead => need to write the Thiele difference equations as a pair of
+ # recursive equations rather than a single recursive formula...)
+ afterDeath = function(advance = NULL, arrears = NULL, start = 0, ..., m = private$m, mCorrection = private$mCorrection, v = private$v) {
+ # General Note: Since the CF vectors can have an arbitrary number of
+ # dimensions, we cannot directly access them via advance[1,..]. Rather,
+ # we have to construct the `[` function manually as a quoted expression,
+ # inserting the required number of dimensions and then evaluating that
+ # expression. This makes this function a little harder to read, but the
+ # performance should not take a hit and it is implemented in a very
+ # general way.
+
+ cfs = list(advance, arrears)
+ # https://stackoverflow.com/a/16896422/920231
+ deflt = cfs[!unlist(lapply(cfs, is.null))][[1]] * 0
+
+ if (missing(advance) || is.null(advance)) advance = deflt
+ if (missing(arrears) || is.null(arrears)) arrears = deflt
+
+ l = max(unlist(lapply(cfs, function(cf) if(!is.null(dim(cf))) dim(cf)[[1]] else length(cf))))
+
+ # TODO: Make sure all CF tensors have the same number of dimensions
+ dms = if (is.null(dim(advance))) length(advance) else dim(advance);
+
+ # Resulting PV tensor has one timestep more than the CF tensors!
+ dms[1] = dms[1] + 1
+ dmNr = if (is.null(dim(advance))) 1 else length(dim(advance))
+
+ # To be able to access the CF tensors in arbitrary dimensions, we
+ # construct the [..] operator manually by quoting it and then inserting
+ # arguments matching the number of dimensions
+ Qadv = Quote(advance[] )[c(1,2,rep(3, dmNr))];
+ Qarr = Quote(arrears[] )[c(1,2,rep(3, dmNr))];
+ Qres = Quote(res[])[c(1,2,rep(3, dmNr))]; # Access the correct number of dimensions
+ QresAss = Quote(res <- tmp)
+
+ VL = function(quoted, time) {
+ eval(quoted %>% `[<-`(3, time))
+ }
+
+ init = VL(Qadv, 1) * 0;
+
+ # assuming advance and arrears have the same dimensions...
+ p = pad0(private$qx$px, l, value=0);
+ q = pad0(private$qx$qx, l, value=0);
+ # TODO: Pad to length l
+ # advance = pad0(advance, l, value = init);
+ # arrears = pad0(arrears, l, value = init);
+
+ # TODO: Replace loop by better way (using Reduce?)
+ res = array(0, dim = dms)
+
+ # Starting value for the recursion:
+ prev = init;
+ prev.dead = init;
+ QresAss[[2]] = Qres %>% `[<-`(3, dms[[1]])
+ eval(QresAss)
+
+ advcoeff = mCorrection$alpha - mCorrection$beta * (1 - v);
+ arrcoeff = mCorrection$alpha - (mCorrection$beta + 1/m) * (1 - v);
+ for (i in l:(start + 1)) {
+ # The actual recursion:
+ tmp = p[i] * v * prev + q[i] * v * prev.dead;
+ # Assign tmp to the slice res[i, ....]
+ QresAss[[2]] = Qres %>% `[<-`(3, i)
+ eval(QresAss)
+ prev = tmp
+ prev.dead = VL(Qadv, i) * advcoeff[i] + VL(Qarr, i) * arrcoeff[i] + v * prev.dead;
+ }
+ VL(Qres, list(1:l))
}
- }
- res
-}
+ ),
+ private = list(
+ qx = data.frame(Alter = c(), qx = c(), ix = c(), px = c()),
+ m = 1,
+ mCorrection = list(alpha = 1, beta = 0),
+ v = 1
+ )
+);
diff --git a/R/InsuranceParameters.R b/R/InsuranceParameters.R
index 72c8f1916dfa0538c898ddfa89cbf17f73cee98a..cd44b442c0ed63d186f18c7b52daccd73eb2f4fe 100644
--- a/R/InsuranceParameters.R
+++ b/R/InsuranceParameters.R
@@ -42,7 +42,7 @@ setCost = function(costs, type, basis = "SumInsured", frequency = "PolicyPeriod"
#' Initialize a cost matrix with dimensions: {CostType, Basis, Period}, where:
#' \describe{
#' \item{CostType:}{alpha, Zillmer, beta, gamma, gamma_nopremiums, unitcosts}
-#' \item{Basis:}{SumInsured, SumPremiums, GrossPremium, NetPremium, Constant}
+#' \item{Basis:}{SumInsured, SumPremiums, GrossPremium, NetPremium, Benefits, Constant}
#' \item{Period:}{once, PremiumPeriod, PremiumFree, PolicyPeriod, CommissionPeriod}
#' }
#' This cost structure can then be modified for non-standard costs using the [setCost()] function.
@@ -100,8 +100,8 @@ setCost = function(costs, type, basis = "SumInsured", frequency = "PolicyPeriod"
initializeCosts = function(costs, alpha, Zillmer, alpha.commission, beta, gamma, gamma.paidUp, gamma.premiumfree, gamma.contract, gamma.afterdeath, gamma.fullcontract, unitcosts, unitcosts.PolicyPeriod) {
if (missing(costs)) {
dimnm = list(
- type = c("alpha", "Zillmer", "beta", "gamma", "gamma_nopremiums", "unitcosts"),
- basis = c("SumInsured", "SumPremiums", "GrossPremium", "NetPremium", "Constant", "Reserve"),
+ type = c("alpha", "Zillmer", "beta", "gamma", "gamma_nopremiums", "unitcosts"),
+ basis = c("SumInsured", "SumPremiums", "GrossPremium", "NetPremium", "Benefits", "Constant", "Reserve"),
frequency = c("once", "PremiumPeriod", "PremiumFree", "PolicyPeriod", "AfterDeath", "FullContract", "CommissionPeriod")
);
costs = array(
diff --git a/R/InsuranceTarif.R b/R/InsuranceTarif.R
index 4bf183ad5ad6374671c3fbbfb6e35abbbc19484e..8021a9cce975d8a3b0660fbdee942cc2a6155139 100644
--- a/R/InsuranceTarif.R
+++ b/R/InsuranceTarif.R
@@ -706,55 +706,34 @@ InsuranceTarif = R6Class(
}
qq = self$getTransitionProbabilities(params, values);
- qx = pad0(qq$qx, values$int$l, value = 1); # After maxAge of the table, use qx=1, also after contract period, just in case
- ix = pad0(qq$ix, values$int$l);
- px = pad0(qq$px, values$int$l);
i = params$ActuarialBases$i;
v = 1/(1 + i);
- benefitFreqCorr = correctionPaymentFrequency(
- i = i, m = params$ContractData$benefitFrequency,
- order = valueOrFunction(params$ActuarialBases$benefitFrequencyOrder, params = params, values = values));
- premiumFreqCorr = correctionPaymentFrequency(
- i = i, m = params$ContractData$premiumFrequency,
- order = valueOrFunction(params$ActuarialBases$premiumFrequencyOrder, params = params, values = values));
-
- pvRefund = calculatePVDeath(px, qx, values$cashFlows$death_GrossPremium, v = v);
- pvRefundPast = calculatePVDeath(
- px, qx,
- values$cashFlows$death_Refund_past,
- v = v) *
+ benefitFreqCorr = correctionPaymentFrequency(
+ i = i, m = params$ContractData$benefitFrequency,
+ order = valueOrFunction(params$ActuarialBases$benefitFrequencyOrder, params = params, values = values));
+ premiumFreqCorr = correctionPaymentFrequency(
+ i = i, m = params$ContractData$premiumFrequency,
+ order = valueOrFunction(params$ActuarialBases$premiumFrequencyOrder, params = params, values = values));
+
+ pvf = PVfactory$new(qx = qq, m = params$ContractData$benefitFrequency, mCorrection = benefitFreqCorr, v = v);
+
+ pvRefund = pvf$death(values$cashFlows$death_GrossPremium);
+ pvRefundPast = pvf$death(values$cashFlows$death_Refund_past) *
(values$cashFlows[,"death_GrossPremium"] - values$cashFlows[,"premiums_advance"]);
pv = cbind(
- premiums = calculatePVSurvival(
- px, qx,
- values$cashFlows$premiums_advance, values$cashFlows$premiums_arrears,
- m = params$ContractData$premiumFrequency, mCorrection = premiumFreqCorr,
- v = v),
- additional_capital = calculatePVSurvival(px, qx, values$cashFlows$additional_capital, 0, v = v),
- guaranteed = calculatePVGuaranteed(
- values$cashFlows$guaranteed_advance, values$cashFlows$guaranteed_arrears,
- m = params$ContractData$benefitFrequency, mCorrection = benefitFreqCorr,
- v = v),
- survival = calculatePVSurvival(
- px, qx,
- values$cashFlows$survival_advance, values$cashFlows$survival_arrears,
- m = params$ContractData$benefitFrequency, mCorrection = benefitFreqCorr,
- v = v),
- death_SumInsured = calculatePVDeath(
- px, qx,
- values$cashFlows$death_SumInsured,
- v = v),
- disease_SumInsured = calculatePVDisease(
- px, qx, ix,
- values$cashFlows$disease_SumInsured, v = v),
+ premiums = pvf$survival(values$cashFlows$premiums_advance, values$cashFlows$premiums_arrears,
+ m = params$ContractData$premiumFrequency, mCorrection = premiumFreqCorr),
+ additional_capital = pvf$survival(advance = values$cashFlows$additional_capital),
+ guaranteed = pvf$guaranteed(values$cashFlows$guaranteed_advance, values$cashFlows$guaranteed_arrears),
+ survival = pvf$survival(values$cashFlows$survival_advance, values$cashFlows$survival_arrears),
+ death_SumInsured = pvf$death(values$cashFlows$death_SumInsured),
+ disease_SumInsured = pvf$disease(values$cashFlows$disease_SumInsured),
death_GrossPremium = pvRefund,
death_Refund_past = pvRefundPast,
death_Refund_future = pvRefund - pvRefundPast,
- death_PremiumFree = calculatePVDeath(
- px, qx,
- values$cashFlows$death_PremiumFree, v = v)
+ death_PremiumFree = pvf$death(values$cashFlows$death_PremiumFree)
);
rownames(pv) <- pad0(rownames(qq), values$int$l);
@@ -772,10 +751,49 @@ InsuranceTarif = R6Class(
}
len = values$int$l;
q = self$getTransitionProbabilities(params, values);
- qx = pad0(q$qx, len);
- px = pad0(q$px, len);
- v = 1/(1 + params$ActuarialBases$i)
- pvc = calculatePVCosts(px, qx, values$cashFlowsCosts, v = v);
+ i = params$ActuarialBases$i
+ v = 1/(1 + i)
+ pvf = PVfactory(qx = q, v = v)
+
+ costs = values$cashFlowsCosts;
+ pvc = costs * 0;
+
+ # survival cash flows (until death)
+ if (any(costs[,,,"survival"] != 0)) {
+ pvc[,,,"survival"] = pvf$survival(advance = costs[,,,"survival"]);
+ }
+ # guaranteed cash flows (even after death)
+ if (any(costs[,,,"guaranteed"] != 0)) {
+ pvc[,,,"guaranteed"] = pvf$guaranteed(advance = costs[,,,"guaranteed"]);
+ }
+
+ # Cash flows only after death
+ if (any(costs[,,,"after.death"] != 0)) {
+ pvc[,,,"after.death"] = pvf$afterDeath(advance = costs[,,,"after.death"]);
+ }
+
+ # Costs based on actual benefits need to be calculated separately
+ # (need to be recalculated with the benefit CFs multiplied by the cost
+ # rate for each type and each year). In most cases, the cost rates will be
+ # constant, but in general we can not assume this => need to re-calculate
+ if (any(values$cashFlowsCosts[,,"Benefits",] != 0)) {
+ bFreq = params$ContractData$benefitFrequency
+ benefitFreqCorr = correctionPaymentFrequency(
+ i = i, m = bFreq,
+ order = valueOrFunction(params$ActuarialBases$benefitFrequencyOrder, params = params, values = values));
+ pvfben = PVfactory$new(qx = q, m = bFreq, mCorrection = benefitFreqCorr, v = v);
+
+ cfCosts = values$cashFlowsCosts[,,"Benefits",];
+ cf = values$cashFlows;
+
+ # Guaranteed + Survival + Death cover + disease
+ pvc[,,"Benefits",] =
+ pvf$guaranteed(cf$guaranteed_advance * cfCosts, cf$guaranteed_arrears * cfCosts) +
+ pvf$survival(cf$survival_advance * cfCosts, cf$survival_arrears * cfCosts) +
+ pvf$death(cf$death_SumInsured * cfCosts) +
+ pvf$disease(cf$disease_SumInsured * cfCosts);
+ }
+
applyHook(hook = params$Hooks$adjustPresentValuesCosts, val = pvc, params = params, values = values, presentValues = presentValues)
},
@@ -1264,8 +1282,9 @@ InsuranceTarif = R6Class(
ZillmerVerteilungCoeff = pad0((0:r)/r, len, 1);
} else {
q = self$getTransitionProbabilities(params, values);
+ pvf = PVfactory$new(qx = q, v = 1/(1 + params$ActuarialBases$i))
# vector of all รค_{x+t, r-t}
- pvAlphaTmp = calculatePVSurvival(q = pad0(q$qx, len), advance = pad0(rep(1,r), len), v = 1/(1 + params$ActuarialBases$i));
+ pvAlphaTmp = pvf$survival(advance = pad0(rep(1,r), len))
ZillmerVerteilungCoeff = (1 - pvAlphaTmp/pvAlphaTmp[[1]]);
}
alphaRefund = ZillmerSoFar - ZillmerVerteilungCoeff * ZillmerTotal;
@@ -1673,8 +1692,15 @@ InsuranceTarif = R6Class(
#' @param ... currently unused
calculatePresentValues = function(cf, params, values) {
len = dim(cf)[1];
- q = self$getTransitionProbabilities(params, values)
- calculatePVSurvival2D(px = pad0(q$px, len), advance = cf, v = 1/(1 + params$ActuarialBases$i));
+ qq = self$getTransitionProbabilities(params, values);
+ i = params$ActuarialBases$i;
+ premFreq = params$ContractData$premiumFrequency;
+ premFreqCorr = correctionPaymentFrequency(
+ i = i, m = premFreq,
+ order = valueOrFunction(params$ActuarialBases$premiumFrequencyOrder, params = params, values = values));
+
+ pvf = PVfactory$new(qx = qq, m = premFreq, mCorrection = premFreqCorr, v = 1/(1 + i));
+ pvf$survival(advance = cf)
},
#' @description Calculate the premium frequency loading, i.e. the surcharge
diff --git a/man/initializeCosts.Rd b/man/initializeCosts.Rd
index ae23aba2dd933c550787b274c369ddbbfea758f8..372e32576ed7eae2f06c93ae3f44c60c0738b870 100644
--- a/man/initializeCosts.Rd
+++ b/man/initializeCosts.Rd
@@ -58,7 +58,7 @@ even if the insured has already dies (for term-fix insurances)}
Initialize a cost matrix with dimensions: {CostType, Basis, Period}, where:
\describe{
\item{CostType:}{alpha, Zillmer, beta, gamma, gamma_nopremiums, unitcosts}
-\item{Basis:}{SumInsured, SumPremiums, GrossPremium, NetPremium, Constant}
+\item{Basis:}{SumInsured, SumPremiums, GrossPremium, NetPremium, Benefits, Constant}
\item{Period:}{once, PremiumPeriod, PremiumFree, PolicyPeriod, CommissionPeriod}
}
This cost structure can then be modified for non-standard costs using the \code{\link[=setCost]{setCost()}} function.
diff --git a/tests/testthat/test-PVfactory.R b/tests/testthat/test-PVfactory.R
new file mode 100644
index 0000000000000000000000000000000000000000..82385867eb7216cefdb7b4dad62d364b7874cdbc
--- /dev/null
+++ b/tests/testthat/test-PVfactory.R
@@ -0,0 +1,90 @@
+test_that("multiplication works", {
+ library(MortalityTables)
+ mortalityTables.load("Austria_Census")
+
+ # Sample data: Austrian population mortality, interest 3%
+ ags = 40:100
+ v = 1/1.03
+ qx = deathProbabilities(mort.AT.census.2011.unisex, ages = ags)
+ qq = data.frame(x = ags, qx = qx, ix = 0, px = 1 - qx)
+ pvf = PVfactory$new(qx = qq, m = 1, v = v)
+ # For invalidity, simply use half the death prob (split qx into qx and ix, leave px unchanged!)
+ qqI = qq
+ qqI$ix = qqI$qx/2
+ qqI$qx = qqI$ix
+ pvfI = PVfactory$new(qx = qqI, m = 1, v = v)
+
+ #################################### #
+ # Cash Flow Definitions
+ #################################### #
+ # Single payment after 10 years
+ cf1 = c(rep(0, 10), 1)
+ # Annual payments for 10 years
+ cfN = rep(1, 10)
+
+ #################################### #
+ # Guaranteed PV
+ #################################### #
+ PV1adv = v^(10:0)
+ PV1arr = v^(11:1)
+ PVNadv = (1 - v^(10:1))/(1-v)
+ PVNarr = v * PVNadv
+
+ # Check basic present values for correctness (one-dimensional CF vector)
+ expect_equal(as.vector(pvf$guaranteed(advance = cf1)), PV1adv)
+ expect_equal(as.vector(pvf$guaranteed(arrears = cf1)), PV1arr)
+ expect_equal(as.vector(pvf$guaranteed(advance = cfN)), PVNadv)
+ expect_equal(as.vector(pvf$guaranteed(arrears = cfN)), PVNarr)
+
+ # Same cash flows, either understood paid in advance at next timestep or in arrears of previous timestep => same present values
+ expect_equal(c(pvf$guaranteed(advance = c(1,cfN))), 1 + c(pvf$guaranteed(arrears = cfN), 0))
+
+ # Check cash flow arrays
+ # PV of single payment is v^(n-t), PV of annuity is (1-v^(n-t))/(1-v)
+ # Use CF array with those two cash flows
+ cf2d = array(c(cf1, cfN, 0), dim = c(length(cf1),2))
+ expect_equal(pvf$guaranteed(advance = cf2d), array(c(PV1adv, PVNadv, 0), dim = c(length(cf1), 2)))
+
+ # two-dimensional cashflows at each time => 3-dimensional tensor
+ cf3d = array(c(cf1, cfN, 0, cf1, cfN, 0, cfN, 0, cf1), dim = c(length(cf1), 2, 3))
+ expect_equal(pvf$guaranteed(advance = cf3d), array(c(PV1adv, PVNadv, 0, PV1adv, PVNadv, 0, PVNadv, 0, PV1adv), dim = c(length(cf1), 2, 3)))
+
+
+
+ #################################### #
+ # Survival PV
+ #################################### #
+ PV1adv.sv = Reduce(`*`, 1-qx[1:10], init = 1, right = TRUE, accumulate = TRUE) * (v^(10:0))
+ PV1arr.sv = head(Reduce(`*`, 1-qx[1:11], init = 1, right = TRUE, accumulate = TRUE), -1) * (v^(11:1))
+ PVNadv.sv = head(Reduce(function(p, pv1) {pv1 * v * p + 1}, 1-qx[1:10], init = 0, right = TRUE, accumulate = TRUE), -1)
+ PVNarr.sv = head(Reduce(function(p, pv1) { v * p * (1 + pv1)}, 1-qx[1:10], init = 0, right = TRUE, accumulate = TRUE), -1)
+
+ # check basic PV
+ expect_equal(as.vector(pvf$survival(advance = cf1)), PV1adv.sv)
+ expect_equal(as.vector(pvf$survival(arrears = cf1)), PV1arr.sv)
+ expect_equal(as.vector(pvf$survival(advance = cfN)), PVNadv.sv)
+ expect_equal(as.vector(pvf$survival(arrears = cfN)), PVNarr.sv)
+
+ # Check cash flow arrays
+ expect_equal(pvf$survival(advance = cf2d), array(c(PV1adv.sv, PVNadv.sv, 0), dim = c(length(cf1), 2)))
+ expect_equal(pvf$survival(arrears = cf2d), array(c(PV1arr.sv, PVNarr.sv, 0), dim = c(length(cf1), 2)))
+
+ # two-dimensional cashflows at each time => 3-dimensional tensor
+ expect_equal(pvf$survival(advance = cf3d), array(c(PV1adv.sv, PVNadv.sv, 0, PV1adv.sv, PVNadv.sv, 0, PVNadv.sv, 0, PV1adv.sv), dim = c(length(cf1), 2, 3)))
+ expect_equal(pvf$survival(arrears = cf3d), array(c(PV1arr.sv, PVNarr.sv, 0, PV1arr.sv, PVNarr.sv, 0, PVNarr.sv, 0, PV1arr.sv), dim = c(length(cf1), 2, 3)))
+
+
+
+ #################################### #
+ # Death PV
+ #################################### #
+ PVN.death = head(Reduce(function(q, pv1) {q * v + (1 - q) * v * pv1}, qx[1:10], init = 0, right = TRUE, accumulate = TRUE), -1)
+ expect_equal(as.vector(pvf$death(benefits = cfN)), PVN.death)
+
+ #################################### #
+ # Disease PV
+ #################################### #
+ # Death and disease probabilities are equal, so the PV should be equal. If death() is implemented correctly, this can detect errors in
+ expect_equal(pvfI$disease(benefits = cfN), pvfI$death(benefits = cfN))
+
+})