For a forthcoming blog post.

simple_mortgage.py

from pylab import *
from scipy.stats import uniform

gamma = 0.01/12.0
t = arange(0,360)

def compute_payout(payout):
    return (payout * exp( - gamma * t)).sum()

def generate_payout():
    payout = ones(shape=(1,360)) * 100
    for k in range(360):
        if uniform().rvs() < 0.01/12.0:
            payout[:,k:] = 0.0
    return payout


N = 100000
values = zeros(shape=(N,), dtype=float)
for i in range(N):
    values[i] = compute_payout(generate_payout())

hist(values, bins=50, normed=True)
xlabel("Total payout")
ylabel("Approximate probability density")
show()