Python Monte Carlo graphs

November 16, 2011 – 10:56 pm

I wanted to plot the paths of my Monte Carlo simulation of Black Scholes in my previous post.

Using iPython started with the –pylab flag, it’s pretty easy.

from math import exp 
from random import gauss
for trials in range(trials):
    for day in range(time):

The code above gives me this nice graph of the stock price “walks” over time:

Things get a little crazier when I make the volatility 3.4 instead of .34:

Using my one line blackscholes function from the previous post, how does the price of the option change as we increase the volatility from .34 to 3.4?

How does the price of the option change as we increase the number of days to expiry?

Black Scholes Monte Carlo in Python

November 16, 2011 – 4:00 pm

The Black Scholes formula is a partial differential equation that can be used to price the present value of an option under certain assumptions. As with everything, you can read all about it in Wikipedia.

It turns that that figuring out the present value of an option is also a good candidate for a Monte Carlo approach. I went through a very simple example of Monte Carlo in my previous post.

Let’s say I have a put option on a particular stock which will pay out at a certain time in the future (the expiry date for European style options, which this approach is restricted to). The payout will be how much less than a given strike price the stock is at the expiry date. If the stock is over that strike price at the expiry date, it’s worth zero.

In the Black Scholes model, stock prices are assumed to fluctuate, or evolve, according to the formula below:

By a process of discretization, we can chop the time interval above into smaller chunks (say 1 day). We can then assume that the brownian motion over this short period will be a normal (gaussian) distribution with a mean of 0 and a variance of the time interval. This is all nicely explained in Wikipedia here and yields the discretized formula:

So turning that into some python code, to find tomorrow’s stock price from today’s price and a given volatility and riskless rate:

>>> from math import exp
>>> from random import gauss
>>> St0=41.75
>>> rate=.0535
>>> vol=.34
>>> delta=exp((rate-.5*vol**2)*(1.0/365)+vol*(1.0/365)**.5*gauss(0,1))
>>> St0+delta

If we wanted to “walk” this price forward for multiple days (say there’s 61 days to expiry), we could add the following to the above:

>>> time=61
>>> delta=lambda x: \
... x*exp((rate-.5*vol**2)*(1.0/365)+vol*(1.0/365)**.5*gauss(0,1))-x
>>> ST=reduce(lambda x,y: x+delta(x),range(time),St0)
>>> ST

So what’s the value of this put option if the strike price is $45?

>>> strike=45
>>> ST=lambda: reduce(lambda x,y: x+delta(x),range(time),St0)
>>> max(0,strike-ST())
>>> max(0,strike-ST())
>>> max(0,strike-ST())
>>> max(0,strike-ST())
>>> max(0,strike-ST())
>>> max(0,strike-ST())

You see that the price of the option varies wildly each time I run this. This is where the Monte Carlo method comes in. I want to run multiple trials of the option value, then average them to get a price that will converge on the Black Scholes price.

A little more python code, and you have it:

>>> trials=10000
>>> reduce(lambda x,y: x+max(0,strike-ST()),range(trials))/trials

So this put option is worth a little over $4 today.

Sorry about the overuse of reduce and lambda functions – I’m kinda smitten by them today for some reason.

Here’s the above nicely formatted:

from math import exp
from random import gauss
delta=lambda x:x*exp((rate-.5*vol**2)*(1.0/365)+vol*(1.0/365)**.5*gauss(0,1))-x
ST=lambda: reduce(lambda x,y: x+delta(x),range(time),St0)
print reduce(lambda x,y: x+max(0,strike-ST()),range(trials))/trials

A friend of mine – Ben – gave me the idea of a function with a custom payout so you can call it with a call, put or other kind of payout (it defaults to a call payout). He’d also I think done a version of it in one line, so here’s my attempt:

from random import gauss
from math import exp
def blackscholes(St0,strike,riskless,vol,time,trials,payout=lambda strike,ST: max(ST-strike,0)):
        return reduce(lambda x,y: x+payout(strike,reduce(lambda x,y: x*exp((riskless-.5*vol**2)*(1.0/365) + vol*(1.0/365)**.5*gauss(0,1)),range(time),St0)),range(trials))/trials

To run this function with the same put example above:

>>> blackscholes(41.75,45,.0535,.34,61,10000,lambda strike,ST:max(strike-ST,0))
>>> 4.195682577752359

By the way, the example above was taken from the Basic Black Scholes book, chapter 5, by Timothy Crack, which I would recommend.

Python Monte Carlo

November 13, 2011 – 9:51 pm

Monte Carlo is a simulation method that can be useful in solving problems that are difficult to solve analytically.

Here’s an interesting application of the technique to estimate the value of pi.

Consider a circular dartboard placed against a square backing, with the sides of the square perfectly touching the circle tangentially at the top, bottom, left and right.

If you throw a dart blindly inside of the square, what are the chances the dart will actually land inside the circle?

Assuming the radius of the circle is r, it is:

= area of circle / area of square

= (pi * r ^ 2) / (4 * r ^ 2)

= pi / 4

So the chances of hitting the circle are pi / 4.   In other words:

pi = 4 * the chances of hitting the circle.

So let’s set up a Monte Carlo simulation to work out the probability of getting the dart inside of the circle.  

To do this, let’s think of the quarter circle below:

Note that the chances of hitting the quarter circle inside of its square is the same as hitting the larger circle inside of the bigger square.   How can we easily simulate hitting the quarter circle?  We can obtain two random numbers between  0 and 1, representing the x and y coordinates of a point with the origin at the circle’s center.  If the distance between this point and the origin is less than 1, then we’re inside the quarter circle.  Otherwise we’re outside it.  If we simulate this enough times, and multiply the percentage of times we’re inside the circle by 4, we should converge on the value of pi.

Here’s a short python script that simulates this:

>>> from random import random
>>> reduce(lambda x,y:x+y,map(lambda x:(random()**2+random()**2)**.5<1,range(100

EDIT 12/4/2011
Thanks curious reader, I guess you’re right, a filter is a little simpler:

>>> from random import random
>>> len(filter(lambda x:(random()**2+random()**2)**.5<1,range(1000000)))/1000000.0*4

Not to be picky, but don’t forget the “.0” at the end of the last 1000000 in your list comprehension solution.

Prime numbers and Numpy – Python

August 27, 2010 – 4:27 pm

I’ve been looking at the element-wise operations in Numpy arrays. The code below generates 100K primes in about 1.3 secs on my laptop.

import numpy
import math
def prime(upto=100):
    return filter(lambda num: (num % numpy.arange(2,1+int(math.sqrt(num)))).all(), range(2,upto+1))
>>> prime()
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

Similiar code that uses a for loop runs about 5 times longer for 100K primes, about 7 secs:

import numpy
import math
def prime2(upto=100):
    return filter(lambda num: numpy.array([num % factor for factor in range(2,1+int(math.sqrt(num)))]).all(), range(2,upto+1))

By not including even numbers, the original code can be optimized to run in half the time again, about .7 secs for 100K primes.

import numpy
import math
def prime3(upto=100):
    return [2]+filter(lambda num: (num % numpy.arange(3,1+int(math.sqrt(num)),2)).all(), range(3,upto+1,2))

The above code is not the fastest implementation, but illustrative of using the element-wise operations in Numpy arrays rather than looping through a list.

We can write faster code by breaking after finding the first modulo equal zero condition. The code below runs in about .4 seconds for 100K primes:

import numpy
import math
def prime4(upto=100):
    for num in range(3,upto+1,2):
        for factor in range(3,1+int(math.sqrt(num)),2):
            if not num % factor: isprime=False; break
        if isprime: primes.append(num)    
    return primes

Ok Scruffy Pete has laid down the challenge in the comments section so I need to get serious about optimization. The code below uses a few Numpy tricks for a zippy implementation and generates 1 million primes in about .03 secs, about 6X faster than Scruffy Pete’s.

import math
import numpy
def prime5(upto):
    for factor in primes[:int(math.sqrt(upto))]:
        if isprime[factor-2]: isprime[factor*2-2::factor]=0
    return primes[isprime]

Scruffy Pete… Ball’s in your court matey!

UPDATE 9/3 – did a half sieve style implementation of the above (see comments) which runs a little faster – under .02 secs for 1 million primes:

import math
import numpy
def prime6(upto):
    for factor in primes[:int(math.sqrt(upto))]:
        if isprime[(factor-2)/2]: isprime[(factor*3-2)/2::factor]=0
    return numpy.insert(primes[isprime],0,2)

I posted this on the stack overflow site (see comments section). They run timings which you can see on:

My impl is there listed under nolfonzo_prime6

Birthday Simulations – Python

August 26, 2010 – 11:38 pm

I’ve wanted for a while to write simple python simulations of the birthday puzzles I’d written about previously. Then, recently, I was wondering what the most efficient way to tell if a list has duplicate elements. It suddenly dawned on me that these two topics were related:

The following code generates birthday samples for a given number of people, then works out the proportion of these samples that contain at least one shared birthday:

import numpy
def sharedbday1(people,simulations=50000):
    samples = numpy.random.random_integers(0,365, size=(simulations, people))
    hasnoshared=map(lambda x: len(x)==len(set(x)),samples)
    return 1-float(sum(hasnoshared))/simulations

Timing this in IPython on my laptop for 50k simulations takes about .5 secs to run:

In [579]: time sharedbday1(23)
CPU times: user 0.56 s, sys: 0.00 s, total: 0.56 s
Wall time: 0.56 s
Out[580]: 0.50678000000000001

Since we’re using Numpy, there’s a cleverer way to see if there’s a shared birthday in each sample. We can sort the array, then compare it element-wise against a slice of itself offset by one – if we find any equivalent elements we’ve found a sample with a shared birthday. With Numpy arrays this is achieved very neatly as operations in Numpy are performed element-wise by default:

import numpy
def sharedbday2(people, simulations=50000):
    samples = numpy.random.random_integers(0,365, size=(simulations, people))
    hasshared  = ((sample[1:]==sample[:-1]).any() for sample in samples)
    return float(sum(hasshared))/simulations

Unfortunately this runs at about half the speed of my original code:

In [582]: time sharedbday2(23)
CPU times: user 1.01 s, sys: 0.00 s, total: 1.01 s
Wall time: 1.02 s
Out[583]: 0.50988

I exchanged some emails with Travis Oliphant of Numpy fame, who suggested the following code which is really neat and I think shows off the power of Numpy in a simple and elegant way:

import numpy
def sharedbday3(people, simulations=50000):
    samples = numpy.random.random_integers(0,365, size=(simulations, people))
    hasshared = (samples[:,1:] == samples[:,:-1]).any(axis=-1)
    return float(sum(hasshared))/simulations

This is fast and shows a speed-up of least 5X on my original code – with a wall time of .09s:

In [585]: time sharedbday3(23)
CPU times: user 0.10 s, sys: 0.00 s, total: 0.10 s
Wall time: 0.09 s
Out[586]: 0.50119999999999998

As Travis explains it:
You create a 2-d array up-front with the simulations along one axis (i.e. every row of the 2-d array is a separate simulation). Then you can use the sample approach used above, but this time letting NumPy do the loop over the rows (down in C). The sort and any methods on NumPy arrays take an axis keyword which determines the dimension over which the operation works. The operation is then repeated implicitly over the other dimensions. So, in the 3rd version shared birthday there are no for loops. The penalty paid is memory usage: a 2-d array instead of throw-away 1-d arrays. But, of course, less memory management is also a major reason for the speed gain.

For completeness, I made some minor modifications to the original code above to simulate the probability that at least one person in a group shares your birthday (rather than of finding a shared birthday between them).

import numpy
def sharemybday(people,simulations=50000):
    samples = numpy.random.random_integers(0,365, size=(simulations, people))
    hasshared=map(lambda x: mybday in x,samples)
    return float(sum(hasshared))/simulations
In [1509]: sharemybday(253)
Out[1509]: 0.50431999999999999

Shortest Path – Python

August 24, 2010 – 2:26 pm

Problem: Write a python function to calculate the shortest path given a start node (or vertex), an end node and a graph. Use Dijkstra’s algorithm. Return the distance from the start node to the end node, as well as the path taken to get there.

The implementation below sticks pretty closely to the algorithm description in the wikipedia entry, which I turned into something a little more pseudo-code-ish to help me implement it:

Initial state:

  1. give nodes two properties – node.visited and node.distance
  2. set node.distance = infinity for all nodes except the start node which is set to zero.
  3. set node.visited = false for all nodes
  4. set current node = start node.

Current node loop:

  1. if current node = end node, finish and return current.distance & path
  2. for all unvisited neighbors, calc their tentative distance (current.distance + edge to neighbor).
  3. if tentative distance < neighbor's set distance, overwrite it.
  4. set current.isvisited = true.
  5. set current = remaining unvisited node with smallest node.distance

Here’s my implementation – it’s recursive, as suggested by the algorithm description:

import sys
def shortestpath(graph,start,end,visited=[],distances={},predecessors={}):
    """Find the shortest path between start and end nodes in a graph"""
    # detect if it's the first time through, set current distance to zero
    if not visited: distances[start]=0
    if start==end:
        # we've found our end node, now find the path to it, and return
        while end != None:
        return distances[start], path[::-1]
    # process neighbors as per algorithm, keep track of predecessors
    for neighbor in graph[start]:
        if neighbor not in visited:
            neighbordist = distances.get(neighbor,sys.maxint)
            tentativedist = distances[start] + graph[start][neighbor]
            if tentativedist < neighbordist:
                distances[neighbor] = tentativedist
    # neighbors processed, now mark the current node as visited
    # finds the closest unvisited node to the start
    unvisiteds = dict((k, distances.get(k,sys.maxint)) for k in graph if k not in visited)
    closestnode = min(unvisiteds, key=unvisiteds.get)
    # now we can take the closest node and recurse, making it current
    return shortestpath(graph,closestnode,end,visited,distances,predecessors)

if __name__ == "__main__":
    graph = {'a': {'w': 14, 'x': 7, 'y': 9},
            'b': {'w': 9, 'z': 6},
            'w': {'a': 14, 'b': 9, 'y': 2},
            'x': {'a': 7, 'y': 10, 'z': 15},
            'y': {'a': 9, 'w': 2, 'x': 10, 'z': 11},
            'z': {'b': 6, 'x': 15, 'y': 11}}
    print shortestpath(graph,'a','a')
    print shortestpath(graph,'a','b')

        (0, ['a'])
        (20, ['a', 'y', 'w', 'b'])

You’ll see I turned the example in wikipedia into a graph for the test case in the code above.

I found this useful about how to represent and process graphs in python: I think this was written by Guido himself.

There’s a bunch of interesting implementations on the web, including two here. More elegant, but it took me a while to understand them.

Cyclic Linked List – Python function to find the node at the start of the loop

July 18, 2010 – 11:01 pm

I’d blogged previously about writing a python function to find out whether a linked list is acyclic or not – it’s worth a read before tackling this one.

Problem: Given the head node of a cyclic (or circular) linked list, write a function to return the node that is at the start of the loop.

A cyclic linked list would look something like this:


In the case above, A is the head node, and C is the node at the start of the loop.

My Solution

Let’s call x the distance from the head to the start of the loop, and c the number of nodes in the loop.

We kick off a fast and a slow pointer from the head, with the fast pointer going twice as fast as the slow pointer (skipping every other node). By the time the slow pointer gets to the start of the loop, the fast pointer will have a head start of x (the fast pointer is always the same distance from the slow pointer as the slow pointer is from the head).

Given its head start of x, when will the fast pointer catch up to the slow pointer? The fast pointer always needs twice the initial distance from the slow pointer to catch up. Moving forward in the loop, the fast pointer is c-x nodes from the slow pointer, so it needs 2(c-x) nodes to catch up. In that time, the slow pointer will have managed to move c-x nodes away from the start (half what the fast pointer travelled).

The key thing to note here is that both pointers are now x nodes away from the start of the loop looking in a forwards direction. So you’ll see that to find the start of the loop the code below gets the fast and slow pointer to meet in the loop, then walks the slow pointer forward x nodes – the distance from the head to the start.

def startOfCyclicLoop(head):

    >>> e1=Element(1)
    >>> e2=Element(2)
    >>> e3=Element(3)
    >>> e4=Element(4)
    >>> e5=Element(5)
    >>> print startOfCyclicLoop(e1)
    >>> print startOfCyclicLoop(e1)
    >>> print startOfCyclicLoop(e1)
    >>> print startOfCyclicLoop(e1)
    >>> print startOfCyclicLoop(e1)
    while fast!=None and!=None:
        if (fast==slow):
            while (fast!=slow):
            return fast
    return  None

class Element:
    def __init__(self,data,next=None):
    def __str__(self):
        return str(

if __name__ == "__main__":
    import doctest

Puzzle – Hands on a Clock

June 6, 2010 – 12:26 am

For a given time, what will be the angle between the hour and minute hands on a clock?

This is rather a long winded way of solving this, so please let me know if you have something shorter.

The angle of the minute hand (from 12 O’clock):
= (minutes / 60) * 360

The angle of the hour hand (from 12 O’clock) is the angle of the hour plus an adjustment for the minutes:
= (hours / 12) * 360 + (minutes / 60) * (1 / 12) * 360

So to get the angle between them we subract the two, and get:
= (hours / 12) * 360 + (minutes / 60) * (1 / 12) * 360 – (minutes / 60) * 360
= (hours * 30) + minutes * .5 – minutes * 6
= hours * 30 – minutes * 5.5

Any angle over 180 we can subtract from 360 to get the shortest angle between them.

Anagrams – Python

June 4, 2010 – 11:56 pm

Write a Python function to find all the anagrams for each of a supplied list of words. For the word dictionary to find the anagrams you could use the one found on Unix or OS X boxes at /usr/share/dict/words. So for example:

>>> for group in anagrams(['dog','cat','horse']):
...     print group
['dog', 'god']
['act', 'cat']
['horse', 'shoer', 'shore']

When called without a parameter, make the function return a list of all the anagram groups in the word dictionary, ordered by the group of anagrams with the most words:

>>> for group in anagrams():
...     if len(group)>6: print len(group), group
9 ['ester', 'estre', 'reest', 'reset', 'steer', 'stere', 'stree', 'terse', 'tsere']
9 ['angor', 'argon', 'goran', 'grano', 'groan', 'nagor', 'orang', 'organ', 'rogan']
9 ['caret', 'carte', 'cater', 'crate', 'creat', 'creta', 'react', 'recta', 'trace']
8 ['leapt', 'palet', 'patel', 'pelta', 'petal', 'plate', 'pleat', 'tepal']
8 ['laster', 'lastre', 'rastle', 'relast', 'resalt', 'salter', 'slater', 'stelar']
7 ['dater', 'derat', 'detar', 'drate', 'rated', 'trade', 'tread']
7 ['easting', 'gainset', 'genista', 'ingesta', 'seating', 'signate', 'teasing']
7 ['darter', 'dartre', 'redart', 'retard', 'retrad', 'tarred', 'trader']
7 ['least',

'setal', 'slate', 'stale', 'steal', 'stela', 'tales']
7 ['atle', 'laet', 'late', 'leat', 'tael', 'tale', 'teal']
7 ['alem', 'alme', 'lame', 'leam', 'male', 'meal', 'mela']
7 ['lapse', 'salep', 'saple', 'sepal', 'slape', 'spale', 'speal']
7 ['armet', 'mater', 'metra', 'ramet', 'tamer', 'terma', 'trame']
7 ['argel', 'ergal', 'garle', 'glare', 'lager', 'large', 'regal']
7 ['aldern', 'darnel', 'enlard', 'lander', 'lenard', 'randle', 'reland']
7 ['alert', 'alter', 'artel', 'later', 'ratel', 'taler', 'telar']
7 ['arist', 'astir', 'sitar', 'stair', 'stria', 'tarsi', 'tisar']
7 ['aril', 'lair', 'lari', 'liar', 'lira', 'rail', 'rial']
7 ['asteer', 'easter', 'reseat', 'saeter', 'seater', 'staree', 'teaser']
7 ['arpent', 'enrapt', 'entrap', 'panter', 'parent', 'pretan', 'trepan']

Can you think of an elegant (dare I say pythonic?) way to implement the anagrams function?

Here’s my attempt:

def anagrams(wordsin=None):
    for word in f.readlines(): 
    if wordsin!=None: 
        return [ana.get(''.join(sorted(wordin)),[]) for wordin in wordsin]   
        return sorted(ana.values(), key=lambda(v): (len(v)),reverse=True)

Puzzle – Birthday pairings

April 30, 2010 – 2:20 pm

Here’s two classic birthday puzzles:

How many people do you need in a room to have a better than 50% chance of finding at least one shared birthday among them?

Like a lot of these problems, it’s easier to first find the chances of the proposition not happening – that is of not finding a shared birthday.

For this problem, I like to imagine that I have a bunch of 365-card decks. I deal one card from each deck. What are the chances that I haven’t dealt the same card after after dealing from r decks? You can see that this is a good analogy for our problem. Looking back at the previous post on dealing a perfect hand may be a good primer for what follows.

The answer to our card analogy is to take the number of ways I can deal the cards (one from each deck) and get no duplicates divided by the total number of ways I can deal the cards. From the first deck, I have 365 ways to deal one card. From the second deck, I have 364 ways to deal one card and not have duplicates. From the third deck, there are 363 ways, and so on until I get to 365-r+1 ways. So the number of ways I can deal cards without duplicates is:

= 365 * 364 * … * (365-r+1)
= 365! / (365-r)!

Now lets work out how many ways there are to deal the cards without restrictions. The answer is simply:

= 365 ^ r

So now, the probability of not dealing any duplicate cards (and also the chances of not finding anyone with a shared birthday) is:

= (365! / (365-r)!) / 365 ^ r

I was a little stumped at how to solve this, except through some trial and error. Even then, 365! seems to make calculators, or even trusty google, blow up. We can use Sterling’s approximation to find that:

ln(365!) = 1792.33142
ln((365-22)!) = 1663.17935
ln((365-23)!) more info

= 1657.34162

So for r=22, our probability of no shared birthdays:
= e ^ (1 792.33142 – 1 663.17935) / 365 ^ 22 = 0.524310203

For r=23:
= e ^ (1792.33142 – 1657.34162) / 365 ^ 23 = 0.492707724

Now remember that this is for finding no shared birthdays. To have at least 1 shared birthday, we need to subract from 1.

So for r=22, the chances of finding at least 1 shared birthday are:
= 1 – 0.524310203 = 0.475689797

For r=23, the chances are:
= 1 – 0.492707724 = 0.507292276

So you can see that it takes at least 23 people in a room to have a greater than 50/50 chance of finding a shared birthday.

By the way, to get the Stirling approximation of ln(365!) I typed the following into google:

365 * ln(365) – (365) + ln(365)/2 + ln (2*pi) / 2)

How many people do you need in a room with you to have a better than 50% chance of finding someone who shares your birthday?

This one’s a little simpler than the first. Again, it’s easier to first find the probability of no one having your birthday. Let’s imagine people walking in to a room where you are. The first person has a 364/365 chance of not having your birthday. And it’s the same for all subsequent people who walk into the room. So the chances of no one sharing your birthday after r people have walked in are:

= (364/365)^r

So the chance that at least one person shares your birthday is:

= 1 – (364/365) ^ r

Let’s find r when the equation is 1/2, representing the at least 50% chance we were looking for:

1/2 = 1 – (364/365) ^ r
(364/365) ^ r = 1/2
r * log (364/365) = log (1/2)
r = log (1/2) / log (364/365)
r = 252.651989

So you need to have at least 253 people with you in a room to have a better than 50/50 chance of finding someone sharing your birthday. It’s quite amazing how different this answer is to the answer for the first part.