diff --git a/R/InsuranceContract.R b/R/InsuranceContract.R
index 0a113c2b9068476cb75ee00b5ac2605e2d471a14..b2717b944244eb7838c1e7572a58009de67e277a 100644
--- a/R/InsuranceContract.R
+++ b/R/InsuranceContract.R
@@ -340,7 +340,7 @@ InsuranceContract = R6Class(
 
     # Add total benefits present value to the PV array. This can only be done after premium calculation, because e.g. premium refund depends on gross premium!
     calculatePresentValuesAllBenefits = function() {
-      pvAllBenefits = self$tarif$presentValueBenefits(presentValues = self$presentValues, premiums = self$premiums, sumInsured = self$sumInsured );
+      pvAllBenefits = self$tarif$presentValueBenefits(presentValues = self$presentValues, presentValuesCosts = self$presentValuesCosts, premiums = self$premiums, sumInsured = self$sumInsured, premiumSum = self$premiumSum );
       self$presentValues = cbind(self$presentValues, pvAllBenefits)
       self$presentValues
     },
@@ -354,7 +354,7 @@ InsuranceContract = R6Class(
     },
 
     calculateReserves = function() {
-      self$reserves = self$tarif$reserveCalculation(premiums=self$premiums, pvBenefits=self$presentValues, pvCosts=self$presentValuesCosts, sumInsured=self$sumInsured, loadings = self$loadings);
+      self$reserves = self$tarif$reserveCalculation(premiums=self$premiums, pvBenefits=self$presentValues, pvCosts=self$presentValuesCosts, sumInsured=self$sumInsured, premiumSum = self$premiumSum, loadings = self$loadings);
     },
 
     premiumAnalysis = function() {
diff --git a/R/InsuranceTarif.R b/R/InsuranceTarif.R
index e6bbbe9822baafe686fce0451f83311f2ee008f7..a3f5b8989cd552f98af70d981d13539aff66003c 100644
--- a/R/InsuranceTarif.R
+++ b/R/InsuranceTarif.R
@@ -248,12 +248,24 @@ InsuranceTarif = R6Class(
       dimnames(res) = list(nm[[1]], colnames)
       res
     },
-    presentValueBenefits = function(presentValues, premiums, sumInsured=1) {
+    presentValueBenefits = function(presentValues, presentValuesCosts, premiums, sumInsured=1, premiumSum=0) {
+      refundAddon = self$premiumRefundLoading;
+      # TODO: Here we don't use the securityLoading parameter => Shall it be used or are these values to be understood without additional security loading?
       benefits.unit    = presentValues[,"survival"] + presentValues[,"death_SumInsured"];
       benefits         = benefits.unit * sumInsured;
-      allBenefits.unit = presentValues[,"survival"] + presentValues[,"death_SumInsured"] + presentValues[,"death_GrossPremium"] * premiums[["unit.gross"]];
+      allBenefits.unit = presentValues[,"survival"] + presentValues[,"death_SumInsured"] + presentValues[,"death_GrossPremium"] * premiums[["unit.gross"]] * (1+refundAddon);
       allBenefits      = allBenefits.unit * sumInsured;
-      cbind(benefits.unit=benefits.unit, benefits=benefits, allBenefits.unit=allBenefits.unit, allBenefits=allBenefits)
+
+      benefitsCosts = presentValuesCosts[,,"SumInsured"]*sumInsured +
+        presentValuesCosts[,,"SumPremiums"] * premiumSum * premiums[["gross"]] +
+        presentValuesCosts[,,"GrossPremium"] * premiums[["gross"]];
+
+      cbind(
+        benefits.unit=benefits.unit,
+        benefits=benefits,
+        allBenefits.unit=allBenefits.unit,
+        allBenefits=allBenefits,
+        benefitsCosts)
     },
 
     getPremiumCoefficients = function(type="gross", coeffBenefits, coeffCosts, ...,
@@ -365,15 +377,21 @@ InsuranceTarif = R6Class(
       list("premiums"=premiums, "coefficients"=coefficients)
     },
 
-    reserveCalculation = function (premiums, pvBenefits, pvCosts, sumInsured=1, ...) {
+    reserveCalculation = function (premiums, pvBenefits, pvCosts, sumInsured=1, premiumSum=0, ...) {
       resNet = pvBenefits[,"allBenefits"]*(1+self$loadings$security) - premiums[["net"]] * pvBenefits[,"premiums"];
-      resZ = pvBenefits[,"allBenefits"]*(1+self$loadings$security) - premiums[["Zillmer"]] * pvBenefits[,"premiums"];
+      BWLminBWP = pvBenefits[,"allBenefits"]*(1+self$loadings$security) - premiums[["net"]] * pvBenefits[,"premiums"];
+      BWZcorr = pvBenefits["0","Zillmer"]/pvBenefits["0", "premiums"]*pvBenefits[,"premiums"];
+      # resZ = pvBenefits[,"allBenefits"]*(1+self$loadings$security) - premiums[["net"]] * pvBenefits[,"premiums"] - pvCosts["0", "Zillmer", "SumInsured"]*premiumSum/pvBenefits["0", "premiums"]*pvBenefits[,"premiums"];
+      resZ=BWLminBWP - BWZcorr;
+
+        #premiums[["Zillmer"]] * pvBenefits[,"premiums"];
       res.gamma    = (pvCosts[,"gamma", "SumInsured"] - pvCosts["0", "gamma", "SumInsured"]/pvBenefits["0", "premiums"]*pvBenefits[,"premiums"])*sumInsured;
+      res.gamma.alt = pvBenefits[,"gamma"] - pvBenefits["0", "gamma"]/pvBenefits["0", "premiums"]*pvBenefits[,"premiums"]
 
       # res.premiumfree =
       # res.gamma.premiumfree =
 
-      res = cbind("net"=resNet, "Zillmer"=resZ, "gamma"=res.gamma#, "Reserve.premiumfree"=res.premiumfree, "Reserve.gamma.premiumfree"=res.gamma.premiumfree);
+      res = cbind("net"=resNet, "BWLminBWP"=BWLminBWP, "BWZcorr"= BWZcorr, "Zillmer"=resZ, "gamma"=res.gamma, "gamma.alt"=res.gamma.alt#, "Reserve.premiumfree"=res.premiumfree, "Reserve.gamma.premiumfree"=res.gamma.premiumfree);
       );
       rownames(res) <- rownames(pvBenefits);
       res