portfolio.py

from pylab import *
from scipy.stats import uniform

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

num_mortgages=50

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

def generate_payout():
    payout = ones(shape=(num_mortgages,360)) * 100
    for i in range(num_mortgages):
        rv = (uniform().rvs((360,)) < 0.01/12.0)
        rv = rv.cumsum()
        payout[i,where(rv > 0)] = 0.0
    return payout


N = 1000
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()