Calculate the probability of all possible cenarios.

policy_probability_list_comprehension.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

import itertools

def benefit(policy_years):
    """Calculate the benefit for the given number of years."""
    
    return sum([policy_years[x] * policy_holders[x] for x in range(len(policy_years))])

def benefit_probability(amount, policy_years):
    """Calculate the probability for given years, if it is
    higher than the amount."""

    bf = benefit(policy_years)

    benefits = [payments_benefit[x] for x in policy_years]
    product = 1
    for bfs in benefits:
        product *= bfs 

    probability = 0

    if bf > amount:
        probability += product

    return probability

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

    probability = 0

    for cenario in itertools.product(range(11), repeat=6):
        probability += benefit_probability(amount, cenario)

    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)))