Retirement Simulation

import locale
import random
import numpy as np

locale.setlocale(locale.LC_ALL, '')
simulations = 1_000_000
stocks_mix = 0.4
withdrawl_rate = 0.04
max_monthly_withdrawl = 20_000
years = 40
recession_at_first_year = 0.5

results = np.zeros((years, simulations))
for s in range(simulations):
    initial = 6_000_000 * recession_at_first_year
    for y in range(years):

        # rates based on last 40 years
        stocks_return = random.normalvariate(0.101, 0.161) 
        bonds_return = random.normalvariate(0.054, 0.031) 
        inflation = random.normalvariate(0.029, 0.015) 
        return_rate = (stocks_mix * stocks_return) + ((1 - stocks_mix) * bonds_return) - inflation

        returns = 0
        withdrawl = 0
        for m in range(12):
            withdrawl += min(initial * (withdrawl_rate / 12), max_monthly_withdrawl)
            initial -= withdrawl
            returns += initial * (return_rate / 12)
            initial += returns
        results[y][s] = withdrawl

for y in range(years):
    dollars = locale.currency(np.median(results[y]), grouping=True)
    print(f'{y}: {dollars}')