Kolmogorov

Kolmogorov using Recursion

This method solves the probability, but requires more computing power than the subsecuent methods because it does not store subsolutions and must recompute them each four times. O(4^n)

#2P10^00 = 1P10^00 + h sum(1P10^01 *mu11^10 + 1P10^02 *mu11^20 - 
#1P10^00 * (mu11^01+mu11^02))
#1P10^02 = 0P10^02 + h sum(0P10^00 *mu11^02 + 1P10^01 *mu11^12 - 
#1P10^02 * (mu11^21+mu11^20))

thpxij <- function(i, j, t, x, h, jMax = 2) {
  sum <- 0
  muSum <- 0
  for (k in 0:jMax) {
    if (k != j) {
      sum <- sum + tpxij(i, k, t, x, h) * mu(k, j, x + t)
      muSum <- muSum + mu(i, k, x + t)
    }
  }
  return(sum - muSum * tpxij(i, j, t, x, h))
}

tpxij <- function(i, j, t, x, h = 1) {
  #print(t)
  if(t <= 0.001) {
    if(i == j) {
       return(1);
     }
    else {
      return(0);
    }
  } else if(i == 2) {
    if(j == 2) {
      return(1)
    } else {
      return(0)
    }
  } else {
    return(tpxij(i, j, t - 1/h, x, h) + 1/h * thpxij(i, j, t - 1/h, x, h))
  }
}

mu <- function(i, j, x) {
  if(i == 0 && j == 1) {
    return(0.0005 * 1.05^x)
  }
  else if (i == 0 && j == 2) {
    return(0.00035*1.075^x)
  } else if (i == 1 && j == 0) {
    return(0.01)
  } else if (i == 1 && j == 2) {
    return(0.0005 + 0.00035*1.075^x)
  } else {
    return(0);
  }
}

Kolmogorov using Memoization (Not debugged)

This method solves top-down, storing a probability once it has been calculated. Unsolved inheritence issues mean that it is not optimized.

tpxijMem <- function(i, j, t, x, h = 1) {
  memo00 <- c(1, rep.int(-1, t*h + 1))
  memo01 <- c(0, rep.int(-1, t*h + 1))
  memo02 <- c(0, rep.int(-1, t*h + 1))
  memo10 <- c(0, rep.int(-1, t*h + 1))
  memo11 <- c(1, rep.int(-1, t*h + 1))
  memo12 <- c(0, rep.int(-1, t*h + 1))
  
  return(tpxMem(i, j, t, x, h = h, memo00, memo01, memo02, memo10, memo11, memo12))
}

tpxMem <- function(i, j, t, x, h = 1, memo00, memo01, memo02, memo10, memo11, memo12) {
  
  if (i == 0 && j == 0 && memo00[t*h + 1] != -1) {
    print("Got here")
    return(memo00[t*h + 1])
  } else if (i == 0 && j == 1 && memo01[t*h + 1] != -1) {
    return(memo01[t*h + 1])
  } else if (i == 0 && j == 2 && memo02[t*h + 1] != -1) {
    return(memo02[t*h + 1])
  } else if (i == 1 && j == 0 && memo10[t*h + 1] != -1) {
    return(memo10[t*h + 1])
  } else if (i == 1 && j == 1 && memo11[t*h + 1] != -1) {
    return(memo11[t*h + 1])
  } else if (i == 1 && j == 2 && memo12[t*h + 1] != -1) {
    return(memo12[t*h + 1])
  }
  
  if(t <= 0.0001) {
    if(i == j) {
       return(1);
     }
    else {
      return(0);
    }
  } else if(i == 2) {
    if(j == 2) {
      return(1)
    } else {
      return(0)
    }
  } else {
    print(memo00, memo01)
    if (i == 0 && j == 0) {
      memo00[t*h + 1] <- tpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12) + 1/h * thpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12)
      return(memo00[t*h + 1])
    } else if (i == 0 && j == 1) {
      memo01[t*h + 1] <- tpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12) + 1/h * thpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12)
      return(memo01[t*h + 1])
    } else if (i == 0 && j == 2) {
      memo02[t*h + 1] <- tpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12) + 1/h * thpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12)
      return(memo02[t*h + 1])
    } else if (i == 1 && j == 0) {
      memo10[t*h + 1] <- tpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12) + 1/h * thpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12)
      return(memo10[t*h + 1])
    } else if (i == 1 && j == 1) {
      memo11[t*h + 1] <- tpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12) + 1/h * thpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12)
      return(memo11[t*h + 1])
    } else if (i == 1 && j == 2) {
      memo12[t*h + 1] <- tpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12) + 1/h * thpxMem(i, j, t - 1/h, x, h, memo00, memo01, memo02, memo10, memo11, memo12)
      return(memo12[t*h + 1])
    }
  }
}


thpxMem <- function(i, j, t, x, h, memo00, memo01, memo02, memo10, memo11, memo12, jMax = 2) {
  sum <- 0
  muSum <- 0
  for (k in 0:jMax) {
    if (k != j) {
      sum <- sum + tpxMem(i, k, t, x, h, memo00, memo01, memo02, memo10, memo11, memo12) * mu(k, j, x + t)
      muSum <- muSum + mu(i, k, x + t)
    }
  }
  return(sum - muSum * tpxMem(i, j, t, x, h, memo00, memo01, memo02, memo10, memo11, memo12))
}

Kolmogorov using Tabulation

This method solves bottom-up, first calculating the probability for one period before solving the next period. No recursion necessary.

tpxijTab <- function(i, j, T, x, h = 1) {
  tab00 <- c(1, rep.int(-1, T*h))
  tab01 <- c(0, rep.int(-1, T*h))
  tab02 <- c(0, rep.int(-1, T*h))
  tab10 <- c(0, rep.int(-1, T*h))
  tab11 <- c(1, rep.int(-1, T*h))
  tab12 <- c(0, rep.int(-1, T*h))
  
  if (i == 2) {
    if (j == 2) {
      return(1)
    } else {
      return(0)
    }
  }
  
  for(t in 2:(T*h + 1)) {
    sum <- 0
    muSum <- 0
    sum <- tab01[t-1] * mu(1, 0, x + (t-2)/h) + tab02[t-1] * mu(2, 0, x + (t-2)/h)
    muSum <- mu(0, 1, x + (t-2)/h) + mu(0, 2, x + (t-2)/h)
    sumMuSum <- sum - muSum * tab00[t-1]
    tab00[t] = tab00[t-1] + 1/h * sumMuSum
    
    sum <- 0
    muSum <- 0
    sum <- tab00[t-1] * mu(0, 1, x + (t-2)/h) + tab02[t-1] * mu(2, 1, x + (t-2)/h)
    muSum <- mu(1, 0, x + (t-2)/h) + mu(1, 2, x + (t-2)/h)
    sumMuSum <- sum - muSum * tab01[t-1]
    tab01[t] = tab01[t-1] + 1/h * sumMuSum
    
    sum <- 0
    muSum <- 0
    sum <- tab00[t-1] * mu(0, 2, x + (t-2)/h) + tab01[t-1] * mu(1, 2, x + (t-2)/h)
    muSum <- mu(2, 0, x + (t-2)/h) + mu(2, 1, x + (t-2)/h)
    sumMuSum <- sum - muSum * tab02[t-1]
    tab02[t] = tab02[t-1] + 1/h * sumMuSum
    
    sum <- 0
    muSum <- 0
    sum <- tab11[t-1] * mu(1, 0, x + (t-2)/h) + tab12[t-1] * mu(2, 0, x + (t-2)/h)
    muSum <- mu(0, 1, x + (t-2)/h) + mu(0, 2, x + (t-2)/h)
    sumMuSum <- sum - muSum * tab10[t-1]
    tab10[t] = tab10[t-1] + 1/h * sumMuSum
    
    sum <- 0
    muSum <- 0
    sum <- tab10[t-1] * mu(0, 1, x + (t-2)/h) + tab12[t-1] * mu(2, 1, x + (t-2)/h)
    muSum <- mu(1, 0, x + (t-2)/h) + mu(1, 2, x + (t-2)/h)
    sumMuSum <- sum - muSum * tab11[t-1]
    tab11[t] = tab11[t-1] + 1/h * sumMuSum
    
    sum <- 0
    muSum <- 0
    sum <- tab10[t-1] * mu(0, 2, x + (t-2)/h) + tab11[t-1] * mu(1, 2, x + (t-2)/h)
    muSum <- mu(2, 1, x + (t-2)/h) + mu(2, 0, x + (t-2)/h)
    sumMuSum <- sum - muSum * tab12[t-1]
    tab12[t] = tab12[t-1] + 1/h * sumMuSum
  }
  
  if (i == 0) {
    if (j == 0) {
      return(tab00[T*h + 1])
    } else if (j == 1) {
      return(tab01[T*h + 1])
    } else if (j == 2) {
      return(tab02[T*h + 1])
    }
  } else if (i == 1) {
    if (j == 0) {
      return(tab10[T*h + 1])
    } else if (j == 1) {
      return(tab11[T*h + 1])
    } else if (j == 2) {
      return(tab12[T*h + 1])
    }
  }
}

Benchmarking

This section calculates how much computing time elapsed.

st <- proc.time()
probability <- tpxijTab(0, 0, 10, 40, h = 12)
et <- proc.time()

out <- paste("Your probability is", probability, "- Calculated in", (et - st)[[3]], "seconds")
print(out)

## [1] "Your probability is 0.873094615638974 - Calculated in 0.156 seconds"

Premium Calculations

If they are healthy at time t, they pay a dollar at the beginning of the month.
If they are disabled at time t, they get 60% of their salary.

premiumValue <- function(T, i = 0, x = 40, h = 1) {
  premium <- 1
  periods <- (1:((T*h)-1))/h
  for (t in periods) {
    premium <- premium + exp(-0.03*t)*tpxijTab(i, i, t, x, h)
  }
  return(premium);
}

benefitValue <- function(T, j = 1, x = 40, h = 1) {
  benefit <- 0
  salary <- 100000
  periods <- (1:(T*h))/h
  for (t in periods) {
    if (t %% 1 == 0) {
      salary <- salary*1.02
    }
    benefit <- benefit + 0.6 * (1/h) * salary * exp(-0.03*t) * tpxijTab(0, j, t, x, h)
  }
  return(benefit)
}

Final answer

40-year old professor buys a DII policy that pays 60% of annual salary while disabled. Salary increases yearly at 2%. Interest is constant at 3%

(ans <- benefitValue(20, h = 12) / premiumValue(20, h = 12))

## [1] 254.6749

Final answer = 254.6749463