Calculate the probability of all possible cenarios.

policy_probability.py

# This problem lies in calculating the probability
# that a certain cenario has to occur. It is
# similar to a Monte Carlo method in that
# it creates several cenarios. However, MC method
# creates cenarios based on the distribution of
# already seen cases while this method brute forces
# all possible cenarios.

# Original post on reddit: https://www.reddit.com/r/vba/comments/b0d07g/brand_new_to_vba_struggling_with_probability/

# The Happy Retirement Insurance Company (HRIC)
# has six elderly policyholders who are all of 
# the same age. If a policyholder is still alive
# at the end of a year, that the policyholder will
# receive their annual benefit. If a policyholder
# is not alive at year end then nothing will be
# paid on their behalf. Here are the annual benefits
# for each of the six policyholders:

policy_holders = [
    50000,
    50000,
    60000,
    60000,
    70000,
    80000
]

# The table bellow shows HRIC's assumptions about 
# the number of payments N that will be paid to any
# randomly selected policyholder. The future lifetimes
# of the policyholders are assumed to be independent,
# i.e. the death or survival of any policyholder has no
# effect on other policyholders.

payments_benefit = [
    0.16,
    0.15,
    0.14,
    0.12,
    0.10,
    0.09,
    0.07,
    0.06,
    0.05,
    0.04,
    0.02
]

# For example, if policyholder number 6 lives for
# four years but dies before the fifth year then
# $80,000 would be paid to policyholder 6 at the
# end of each year for four years. The probability
# of this $320,000 payout for policyholder 6 is 0.10.

# Questions:
# 1. What the probability of paying 1,200,000 or more
# 2. =    =   =           =  =      1,400,000 =  =
# 3. =    =   =           =  =      1,700,000 =  =

# Answers:
# 1. 0.545043
# 2. 0.366003
# 3. 0.158900

def benefit(a, b, c, d, e, f):
    """Calculate the benefit for the given number of years."""
    
    result = (a * policy_holders[0] +
                b * policy_holders[1] +
                c * policy_holders[2] +
                d * policy_holders[3] +
                e * policy_holders[4] +
                f * policy_holders[5])

    return result

def benefit_probability(amount, a, b, c, d, e, f):
    """Calculate the probability for given years, if it is
    higher than the amount."""

    bf = benefit(a, b, c, d, e, f)

    if bf > amount:
        probability = (payments_benefit[a] * 
                        payments_benefit[b] *
                        payments_benefit[c] *
                        payments_benefit[d] *
                        payments_benefit[e] *
                        payments_benefit[f])
    else:
        probability = 0
    
    return probability

def calculate_probability(amount):
    """Generate all possible cases for each policyholder
    and the number of years (number of payments) from 
    0 (minimum number) to 11 (maximum number)."""

    probability = 0

    # Policyholder 1 ...
    for a in range(11):
        for b in range(11):
            for c in range(11):
                for d in range(11):
                    for e in range(11):
                        for f in range(11):
                            probability += benefit_probability(amount, a, b, c, d, e, f)

    return probability

if __name__ == "__main__":
    question_1 = calculate_probability(1200000)
    question_2 = calculate_probability(1400000)
    question_3 = calculate_probability(1700000)

    # Compare results to the expected results:
    print("Question 1:\n\t- Expected: {}\n\t- Obtained: {}".format(0.545043, round(question_1, 6)))
    print("Question 1:\n\t- Expected: {}\n\t- Obtained: {}".format(0.366603, round(question_2, 6)))
    print("Question 1:\n\t- Expected: {}\n\t- Obtained: {}".format(0.158900, round(question_3, 6)))