diff --git a/R/InsuranceContract.R b/R/InsuranceContract.R
index 6b3121f83a80607873848c017cf8660a9b5e8794..f9eab7dcf5875ec8425b90f07204f43c809dbe43 100644
--- a/R/InsuranceContract.R
+++ b/R/InsuranceContract.R
@@ -1,16 +1,15 @@
 library(R6)
 library(openxlsx);
-# require(xlsx)
 
-InsuranceContract = R6Class(
-  "InsuranceContract",
+InsuranceContract = R6Class("InsuranceContract",
+
   public = list(
     tarif = NA,
 
     #### Contract settings
     params = list(
       sumInsured = 1,
-      premiumWaiver= 0,
+      premiumWaiver = 0,
       YOB = NA,
       age = NA,
       policyPeriod = Inf,
@@ -24,11 +23,14 @@ InsuranceContract = R6Class(
       premiumFrequency = 1,
       benefitFrequency = 1, # Only for annuities!
 
-      loadings = list()     # Allow overriding the tariff-defined loadings (see the InsuranceTariff class for all possible names)
+      loadings = list(),     # Allow overriding the tariff-defined loadings (see the InsuranceTariff class for all possible names)
+      surrenderPenalty = TRUE,  # Set to FALSE after the surrender penalty has been applied once, e.g. on a premium waiver
+      alphaRefunded = FALSE     # Alpha costs not yet refunded (in case of contract changes)
     ),
 
     #### Caching values for this contract, initialized/calculated when the object is created
     values = list(
+      basicData = NA,
       transitionProbabilities = NA,
 
       cashFlowsBasic = NA,
@@ -50,18 +52,15 @@ InsuranceContract = R6Class(
     ),
 
     #### Keeping the history of all contract changes during its lifetime
-    history = list(
-
-    ),
+    history = list(),
 
 
     #### The code:
 
-    initialize = function(tarif, age, policyPeriod,
-                          premiumPeriod = policyPeriod, sumInsured = 1,
+    initialize = function(tarif, age, sumInsured = 1,
+                          policyPeriod, premiumPeriod = policyPeriod, guaranteed = 0,
                           ...,
                           loadings = list(),
-                          guaranteed = 0,
                           premiumPayments = "in advance", benefitPayments = "in advance",
                           premiumFrequency = 1, benefitFrequency = 1,
                           deferral = 0, YOB = 1975) {
@@ -83,102 +82,128 @@ InsuranceContract = R6Class(
       if (!missing(guaranteed))       self$params$guaranteed = guaranteed;
       if (!missing(loadings))         self$params$loadings = loadings;
 
-      self$recalculate();
+      self$calculateContract();
     },
 
     addHistorySnapshot = function(time=0, comment="Initial contract values", type="Contract", params=self$params, values = self$values) {
-      self$history = c(self$history,
-                       list("time"=time, "comment"=comment, "type"=type, "params"=params, "values"=values));
+      self$history = rbind(self$history,
+                       list(time=list("time"=time, "comment"=comment, "type"=type, "params"=params, "values"=values)));
     },
 
-    recalculate = function() {
-      self$determineTransitionProbabilities();
-      self$determineCashFlows();
-      self$calculatePresentValues();
-      self$calculatePremiums();
-      self$updatePresentValues(); # Update the cash flows and present values with the values of the premium
-
-      self$calculateAbsCashFlows();
-      self$calculateAbsPresentValues();
-      self$calculateReserves();
-      self$premiumAnalysis();
-      self$addHistorySnapshot(0, "Initial contract values", type="Contract", params=self$params, values = self$values);
-    },
-
-
-    determineTransitionProbabilities = function() {
-      self$values$transitionProbabilities = do.call(self$tarif$getTransitionProbabilities, self$params);
-      self$values$transitionProbabilities
-    },
+    calculateContract = function() {
+      self$values$transitionProbabilities = self$determineTransitionProbabilities();
 
-    determineCashFlows = function() {
-      self$values$cashFlowsBasic = do.call(self$tarif$getBasicCashFlows, self$params);
-      self$values$cashFlows = do.call(self$tarif$getCashFlows, c(self$params, self$values));
-      self$values$premiumSum = sum(self$values$cashFlows$premiums_advance + self$values$cashFlows$premiums_arrears);
-      self$values$cashFlowsCosts = do.call(self$tarif$getCashFlowsCosts, c(self$params, self$values));
-      list("benefits"= self$values$cashFlows, "costs"=self$values$cashFlowCosts, "premiumSum" = self$values$premiumSum)
-    },
+      self$values$cashFlowsBasic = self$determineCashFlowsBasic();
+      self$values$cashFlows = self$determineCashFlows();
+      self$values$premiumSum = self$determinePremiumSum();
+      self$values$cashFlowsCosts = self$determineCashFlowsCosts();
 
-    calculatePresentValues = function() {
-str(self$values);
-      self$values$presentValues = do.call(self$tarif$presentValueCashFlows,
-                                   c(self$params, self$values));
-      self$values$presentValuesCosts = do.call(self$tarif$presentValueCashFlowsCosts,
-                                        c(self$params, self$values));
-      list("benefits" = self$values$presentValues, "costs" = self$values$presentValuesCosts)
-    },
+      self$values$presentValues = self$calculatePresentValues();
+      self$values$presentValuesCosts = self$calculatePresentValuesCosts();
 
-    calculatePremiums = function() {
       # the premiumCalculation function returns the premiums AND the cofficients,
       # so we have to extract the coefficients and store them in a separate variable
-      res = do.call(self$tarif$premiumCalculation, c(self$params, self$values));
+      res = self$calculatePremiums();
       self$values$premiumCoefficients = res[["coefficients"]];
       self$values$premiums = res[["premiums"]]
-      self$values$premiums
-    },
 
-    updatePresentValues = function() {
-      pvAllBenefits = do.call(self$tarif$presentValueBenefits, c(self$params, self$values));
+      # Update the cash flows and present values with the values of the premium
+      pvAllBenefits = self$calculatePresentValuesBenefits()
       self$values$presentValues = cbind(self$values$presentValues, pvAllBenefits)
-      self$values$presentValue
-    },
 
-    calculateAbsCashFlows = function() {
-      self$values$absCashFlows = do.call(self$tarif$getAbsCashFlows, c(self$params, self$values));
-      self$values$absCashFlows
-    },
+      self$values$absCashFlows = self$calculateAbsCashFlows();
+      self$values$absPresentValues = self$calculateAbsPresentValues();
+      self$values$reserves = self$calculateReserves();
+      self$values$basicData = self$getBasicDataTimeseries()
+      self$values$premiumComposition = self$premiumAnalysis();
 
-    calculateAbsPresentValues = function() {
-      self$values$absPresentValues = do.call(self$tarif$getAbsPresentValues, c(self$params, self$values));
-      self$values$absPresentValues
+      self$addHistorySnapshot(0, "Initial contract values", type="Contract", params=self$params, values = self$values);
     },
 
-    calculateReserves = function() {
-      self$values$reserves = do.call(self$tarif$reserveCalculation, c(self$params, self$values));
-      self$values$reserves
+    determineTransitionProbabilities = function(contractModification=NULL) {
+      do.call(self$tarif$getTransitionProbabilities, c(self$params, self$values, list(contractModification=contractModification)));
     },
-
-    premiumAnalysis = function() {
-      self$values$premiumComposition = do.call(self$tarif$premiumDecomposition, c(self$params, self$values));
-      self$values$premiumComposition
+    determineCashFlowsBasic = function(contractModification=NULL) {
+      do.call(self$tarif$getBasicCashFlows, self$params);
+    },
+    determineCashFlows = function(contractModification=NULL) {
+      do.call(self$tarif$getCashFlows, c(self$params, self$values, list(contractModification=contractModification)));
+    },
+    determinePremiumSum = function(contractModification=NULL) {
+      sum(self$values$cashFlows$premiums_advance + self$values$cashFlows$premiums_arrears);
+    },
+    determineCashFlowsCosts = function(contractModification=NULL) {
+      do.call(self$tarif$getCashFlowsCosts, c(self$params, self$values, list(contractModification=contractModification)));
+    },
+    calculatePresentValues = function(contractModification=NULL) {
+      do.call(self$tarif$presentValueCashFlows, c(self$params, self$values, list(contractModification=contractModification)));
+    },
+    calculatePresentValuesCosts = function(contractModification=NULL) {
+      do.call(self$tarif$presentValueCashFlowsCosts, c(self$params, self$values, list(contractModification=contractModification)));
+    },
+    calculatePremiums = function(contractModification=NULL) {
+      do.call(self$tarif$premiumCalculation, c(self$params, self$values, list(contractModification=contractModification)));
+    },
+    calculatePresentValuesBenefits = function(contractModification=NULL) {
+      do.call(self$tarif$presentValueBenefits, c(self$params, self$values, list(contractModification=contractModification)));
+    },
+    calculateAbsCashFlows = function(contractModification=NULL) {
+      do.call(self$tarif$getAbsCashFlows, c(self$params, self$values, list(contractModification=contractModification)));
+    },
+    calculateAbsPresentValues = function(contractModification=NULL) {
+      do.call(self$tarif$getAbsPresentValues, c(self$params, self$values, list(contractModification=contractModification)));
+    },
+    calculateReserves = function(contractModification=NULL) {
+      do.call(self$tarif$reserveCalculation, c(self$params, self$values, list(contractModification=contractModification)));
+    },
+    premiumAnalysis = function(contractModification=NULL) {
+      do.call(self$tarif$premiumDecomposition, c(self$params, self$values, list(contractModification=contractModification)));
+    },
+    getBasicDataTimeseries = function(contractModification=NULL) {
+      do.call(self$tarif$getBasicDataTimeseries, c(self$params, self$values, list(contractModification=contractModification)));
     },
 
     # Premium Waiver: Stop all premium payments at time t
     # the SumInsured is determined from the available
     premiumWaiver = function (t) {
       newSumInsured = self$values$reserves[[toString(t), "PremiumFreeSumInsured"]];
-      self$premiumWaiver = TRUE;
-      self$recalculatePremiumFreeSumInsured(t=t, SumInsured=newSumInsured)
+      self$params$premiumWaiver = TRUE;
+      self$params$surrenderPenalty = FALSE; # Surrencer penalty has already been applied, don't apply a second time
+      self$params$alphaRefunded = TRUE;     # Alpha cost (if applicable) have already been refunded partially, don't refund again
+
+      self$params$sumInsured = newSumInsured;
+
+      self$values$cashFlowsBasic = mergeValues(starting=self$values$cashFlowsBasic, ending=self$determineCashFlowsBasic(t), t=t);
+      self$values$cashFlows = mergeValues(starting=self$values$cashFlows, ending=self$determineCashFlows(t), t=t);
+      # Premium sum is not affected by premium waivers, i.e. everything depending on the premium sum uses the original premium sum!
+      # self$values$premiumSum = self$determinePremiumSum();
+      self$values$cashFlowsCosts = mergeValues3D(starting=self$values$cashFlowsCosts, ending=self$determineCashFlowsCosts(t), t=t);
+
+      pv = self$calculatePresentValues(t);
+      pvc = self$calculatePresentValuesCosts(t);
+      self$values$presentValuesCosts = mergeValues3D(starting=self$values$presentValuesCosts, ending=pvc, t=t);
 
-      # TODO: Update cashflows from t onwards
-      # TODO: Update present values from t onwards
-      # TODO: Update reserves from t onwards
+      # TODO:
+      # the premiumCalculation function returns the premiums AND the cofficients,
+      # so we have to extract the coefficients and store them in a separate variable
+      # res = self$calculatePremiums(t);
+      # self$values$premiumCoefficients = mergeValues(starting=self$values$premiumCoefficients, ending=res[["coefficients"]], t=t);
+      # self$values$premiums = mergeValues(starting= = res[["premiums"]]
 
+      # Update the cash flows and present values with the values of the premium
+      pvAllBenefits = self$calculatePresentValuesBenefits()
+      self$values$presentValues = mergeValues(starting=self$values$presentValues, ending=cbind(pv, pvAllBenefits), t=t);
 
-      self$addHistorySnapshot(t=t, comment=sprintf("Premium waiver at time %d", t), type="PremiumWaiver", params=self$params, values=self$values);
+      self$values$absCashFlows = mergeValues(starting=self$values$absCashFlows, ending=self$calculateAbsCashFlows(t), t=t);
+      self$values$absPresentValues = mergeValues(starting=self$values$absPresentValues, ending=self$calculateAbsPresentValues(t), t=t);
+      self$values$reserves = mergeValues(starting=self$values$reserves, ending=self$calculateReserves(t), t=t);
+      self$values$basicData = mergeValues(starting=self$values$basicData, ending=self$getBasicDataTimeseries(t), t=t);
+      self$values$premiumComposition = mergeValues(starting=self$values$premiumComposition, ending=self$premiumAnalysis(t), t=t);
+
+      self$addHistorySnapshot(time=t, comment=sprintf("Premium waiver at time %d", t), type="PremiumWaiver", params=self$params, values=self$values);
     },
 
     dummy=NULL
   )
 );
-
+# InsuranceContract$debug("premiumWaiver")
diff --git a/R/InsuranceTarif.R b/R/InsuranceTarif.R
index f0a463c60088ddd822f2392531508b0ac68574fe..72eb8087d99658f0f8543b4eefb178e5ca260d14 100644
--- a/R/InsuranceTarif.R
+++ b/R/InsuranceTarif.R
@@ -1,7 +1,6 @@
 library(R6)
-library(lifecontingencies)
+# library(lifecontingencies)
 library(objectProperties)
-library(foreach)
 
 TariffTypeEnum = setSingleEnum("TariffType", levels = c("annuity", "wholelife", "endowment", "pureendowment", "terme-fix", "dread-disease"))
 PaymentTimeEnum = setSingleEnum("PaymentTime", levels = c("in advance", "in arrears"))
@@ -128,7 +127,7 @@ InsuranceTarif = R6Class(
     },
     getBasicCashFlows = function(age, ..., guaranteed = 0, policyPeriod = inf, deferral = 0, maxAge = getOmega(self$mortalityTable)) {
       maxlen = min(maxAge - age, policyPeriod);
-      cf = list(
+      cf = data.frame(
         guaranteed = rep(0, maxlen+1),
         survival = rep(0, maxlen+1),
         death = rep(0, maxlen+1),
@@ -153,7 +152,7 @@ InsuranceTarif = R6Class(
       cf
     },
 
-    getCashFlows = function(age, ..., premiumPayments = "in advance", benefitPayments = "in advance", guaranteed = 0, policyPeriod=Inf, premiumPeriod = policyPeriod, deferral=0, maxAge = getOmega(self$mortalityTable), cashFlowsBasic = NULL) {
+    getCashFlows = function(age, ..., premiumPayments = "in advance", benefitPayments = "in advance", guaranteed = 0, policyPeriod=Inf, premiumPeriod = policyPeriod, deferral=0, maxAge = getOmega(self$mortalityTable), cashFlowsBasic = NULL, premiumWaiver = FALSE) {
       if (missing(cashFlowsBasic)) {
         cashFlowsBasic = self$getBasicCashFlows(age = age, ..., guaranteed = guaranteed,
               policyPeriod = policyPeriod, deferral = deferral, maxAge = maxAge);
@@ -177,11 +176,13 @@ InsuranceTarif = R6Class(
       );
 
       # Premiums:
-      premiums = pad0(rep(1, min(premiumPeriod, policyPeriod)), cflen);
-      if (premiumPayments == "in advance") {
-        cf$premiums_advance = premiums;
-      } else {
-        cf$premiums_arrears = premiums;
+      if (!premiumWaiver) {
+        premiums = pad0(rep(1, min(premiumPeriod, policyPeriod)), cflen);
+        if (premiumPayments == "in advance") {
+          cf$premiums_advance = premiums;
+        } else {
+          cf$premiums_arrears = premiums;
+        }
       }
 
       # Survival Benefits
@@ -210,7 +211,7 @@ InsuranceTarif = R6Class(
       cf
     },
 
-    getCashFlowsCosts = function(age, ..., policyPeriod=Inf, premiumPeriod = policyPeriod, maxAge = getOmega(self$mortalityTable)) {
+    getCashFlowsCosts = function(age, ..., policyPeriod=Inf, premiumPeriod = policyPeriod, premiumWaiver = FALSE, maxAge = getOmega(self$mortalityTable)) {
       maxlen = min(maxAge - age, policyPeriod)+1;
       policyPeriod = min(maxAge - age, policyPeriod);
       premiumPeriod = min(policyPeriod, premiumPeriod);
@@ -230,6 +231,11 @@ InsuranceTarif = R6Class(
       for (i in 1:policyPeriod) {
         cf[i,,] = cf[i,,] + self$costs[,,"PolicyPeriod"];
       }
+
+      # After premiums are waived, use the gamma_nopremiums instead of gamma:
+      if (premiumWaiver) {
+        cf[,"gamma",] = cf[,"gamma_nopremiums",];
+      }
       cf
     },
 
@@ -245,8 +251,6 @@ InsuranceTarif = R6Class(
       pvRefund = calculatePVDeath (px, qx, cashFlows$death_GrossPremium, v=self$v);
       pvRefundPast = calculatePVDeath (px, qx, cashFlows$death_Refund_past, v=self$v) * (cashFlows[,"death_GrossPremium"]-cashFlows[,"premiums_advance"]);
 
-str(px/px);
-str(qx*0);
       pv = cbind(
         premiums = calculatePVSurvival (px, qx, cashFlows$premiums_advance, cashFlows$premiums_arrears, m=premiumFrequency, mCorrection=premiumFrequencyCorrection, v=self$v),
         guaranteed = calculatePVGuaranteed (cashFlows$guaranteed_advance, cashFlows$guaranteed_arrears, m=benefitFrequency, mCorrection=benefitFrequencyCorrection, v=self$v),
@@ -269,7 +273,6 @@ str(qx*0);
       qx = pad0(q$q, len);
       px = pad0(q$p, len);
 
-# str(cashFlowsCosts);
       pvc = calculatePVCosts(px, qx, cashFlowsCosts, v=self$v);
       pvc
     },
@@ -459,7 +462,7 @@ str(qx*0);
       list("premiums"=premiums, "coefficients"=coefficients)
     },
 
-    reserveCalculation = function (premiums, absPresentValues, absCashFlows, sumInsured=1, premiumSum=0, policyPeriod = 1, age = 0, ..., loadings=list()) {
+    reserveCalculation = function (premiums, absPresentValues, absCashFlows, sumInsured=1, premiumSum=0, policyPeriod = 1, age = 0, ..., reserves = c(), loadings=list(), surrenderPenalty = TRUE) {
       # Merge a possibly passed loadings override with the defaults of this class:
       loadings = self$getLoadings(loadings=loadings);
       # Net, Zillmer and Gross reserves
@@ -505,7 +508,10 @@ str(qx*0);
       # The surrender value functions can have arbitrary form, so we store a function
       # here in the tarif and call that, passing the reduction reserve as
       # starting point, but also all reserves, cash flows, premiums and present values
-      if (!is.null(self$surrenderValueCalculation)) {
+      if (!surrenderPenalty) {
+        # No surrender penalty any more (has already been applied to the first contract change!)
+        surrenderValue = resReduction;
+      } else if (!is.null(self$surrenderValueCalculation)) {
         surrenderValue = self$surrenderValueCalculation(
           resReduction, reserves=res, premiums=premiums, absPresentValues=absPresentValues,
           absCashFlows=absCashFlows, sumInsured=sumInsured, premiumSum=premiumSum,
@@ -530,14 +536,14 @@ str(qx*0);
 
     getBasicDataTimeseries = function(premiums, reserves, absCashFlows, absPresentValues, sumInsured=1, policyPeriod, premiumPeriod, ...) {
       res=cbind(
-        "PremiumPayment" = c(rep(1, premiumPeriod), rep(0, policyPeriod-premiumPeriod)),
-        "SumInsured" = rep(sumInsured, policyPeriod),
+        "PremiumPayment" = c(rep(1, premiumPeriod), rep(0, policyPeriod-premiumPeriod+1)),
+        "SumInsured" = c(rep(sumInsured, policyPeriod), 0),
         "Premiums" = absCashFlows$premiums_advance + absCashFlows$premiums_arrears,
-        "InterestRate" = rep(self$i, policyPeriod),
-        "PolicyDuration" = rep(policyPeriod, policyPeriod),
-        "PremiumPeriod" = rep(premiumPeriod, policyPeriod)
+        "InterestRate" = rep(self$i, policyPeriod+1),
+        "PolicyDuration" = rep(policyPeriod, policyPeriod+1),
+        "PremiumPeriod" = rep(premiumPeriod, policyPeriod+1)
       );
-      rownames(res) = 0:(policyPeriod-1);
+      rownames(res) = 0:policyPeriod;
       res
     },