Proper simulation of annuities

simulate_annuity.py

#!/usr/bin/env python

def simulate_annuity(amount, annuity, month_mean_interest, month_scale_interest):
    import numpy as np
    
    months = []
    
    while amount > 0:
        months.append(amount)
        if month_scale_interest:
            new_interest = np.random.normal(month_mean_interest, month_scale_interest, 1)[0]
        else:
            new_interest = month_mean_interest
            
        amount = amount * math.pow(new_interest, 1./12) - annuity
    return months

if __name__ == "__main__":
    import matplotlib.pyplot as plt
    import math
    
    l = []
    
    for i in range(1000):
      r = simulate_annuity(200000., 1000., 1.03, .2)
    
      plt.plot(r)
      l.append(len(r))
    
    plt.figure()
    plt.hist(l)
    plt.show()
    print sum(l)/len(l)