Revised Monte Carlo simulation of the payoffs to angel investing

angel_sim2.py

#!/usr/bin/python
#
# Monte Carlo simulation of the payoffs to angel investing.
#
# Assume a pool of N different investors, each investing in D deals,
# with a fixed distribution of payoffs per deal.  Randomly simulate
# each investor's combined payoff, then compute the mean and std dev
# of all payoffs in the overall pool.
#
# This gives an individual angel an idea of what kind of payoff &
# variance to expect from investing in a certain # of deals.
#
# Written by Jeff Miller @jmillerinc.
#
# See: http://jmillerinc.com/2010/04/28/angel-investing-simulation-part-2
#
# Payoff distribution from http://www.gabrielweinberg.com/blog/2010/03/
# angel-investing-portfolio-scenario-planner-spreadsheet.html

import math
import optparse
import random
import sys
import unittest


DEFAULT_NUM_INVESTORS = 10000
DEFAULT_NUM_DEALS = 20
PROBS   = [.50, .20, .15, .13, .02]
PAYOFFS = [  0,   1,   3,  10,  20]


def random_draw(cum_probs):
    """Randomly draw from a discrete distribution with the given
    cumulative probabilities. Return the index of the probability
    bucket that was drawn."""
    z = random.random()
    for i, p in enumerate(cum_probs):
        if z <= p:
            return i
    raise Exception('Execution should never reach this point: ' +
                    'cumulative probabilities are invalid.')


def cumulative_probabilities(probs):
    """Given a discrete distribution as a list of probabilities,
    return the list of cumulative probabilities."""
    result = [sum(probs[0:i+1]) for i in range(len(probs))]
    return result


class TestCumulativProbabilities(unittest.TestCase):
    def test(self):
        test_cases = (
            ([1], [1]),
            ([.5, .5], [.5, 1]),
            ([0, .1, .2, .3, .4], [0, .1, .3, .6, 1]),
            )
        for input, expected_output in test_cases:
            output = cumulative_probabilities(input)
            assert abs(output[-1] - 1.0) < 1e-6
            self.assertEqual(len(output), len(expected_output))
            for i in range(len(output)):
                assert abs(output[i] - expected_output[i]) < 1e-10


def median(x):
    """Return the median of a list of numbers."""
    n = len(x)
    y = sorted(x)
    if (n % 2) == 0:
        return 0.5 * (y[n/2] + y[n/2-1])
    else:
        return y[n/2]


class TestMedian(unittest.TestCase):
    def test(self):
        test_cases = (
            ([3], 3),
            ([3, 1.2], 2.1),
            ([3, 1.2, 4], 3),
            ([0, 0, 0, 99, 100, 100, 100], 99),
            ([0, 0, 99, 100, 100, 100], 99.5),
            )
        for input, expected_output in test_cases:
            self.assertEqual(median(input), expected_output)


def simulate(num_investors, num_deals):
    """Simulate N investors doing D deals each."""
    cum_probs = cumulative_probabilities(PROBS)
    deal_size = 1.0/float(num_deals)
    payoffs = []
    for i in range(num_investors):
        payoffs.append(sum(deal_size * PAYOFFS[random_draw(cum_probs)] for j in range(num_deals)))

    n = len(payoffs)
    assert n == num_investors
    summ = sum(payoffs)
    summ2 = sum(x*x for x in payoffs)
    mean_payoff = summ/n
    median_payoff = median(payoffs)
    std_payoff = math.sqrt((summ2 - summ*summ/n)/(n-1))
    idx_less_than_10x_payoff = [i for (i, p) in enumerate(cum_probs) if PAYOFFS[i] < 10]
    prob_less_than_10x_payoff = cum_probs[idx_less_than_10x_payoff[-1]]
    prob_no_hit = pow(prob_less_than_10x_payoff, num_deals)
    print 'N = %5d  D = %3d  mean = %8.4f  median = %8.4f  std = %8.4f  nohit = %8.4f' % \
        (num_investors, num_deals, mean_payoff, median_payoff, std_payoff, prob_no_hit)


if __name__ == '__main__':
    option_parser = optparse.OptionParser()
    option_parser.add_option('-n', type='int', dest='num_investors', default=DEFAULT_NUM_INVESTORS,
                             metavar='N', help='Simulate N investors (default %d)' % DEFAULT_NUM_INVESTORS)
    option_parser.add_option('-d', type='int', dest='num_deals', default=DEFAULT_NUM_DEALS,
                             metavar='D', help='Simulate D deals per investor (default %d)' % DEFAULT_NUM_DEALS)
    option_parser.add_option('--batch', dest='batch', action='store_true', default=False,
                             help='Simulate a sequence of deal sizes.')
    option_parser.add_option('--unittest', dest='unittest', action='store_true', default=False,
                             help='Run unit tests and exit.')
    options, args = option_parser.parse_args()

    if options.unittest:
        unittest.main(argv=sys.argv[0:1])
        sys.exit(0)

    if options.batch:
        num_deals = 1
        while num_deals <= 100:
            simulate(options.num_investors, num_deals)
            if num_deals < 20:
                num_deals += 1
            elif num_deals < 50:
                num_deals += 5
            else:
                num_deals += 10
    else:
        simulate(options.num_investors, options.num_deals)