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 42.747744977337405 >>>
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 38.89640199216139
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()) 6.323401387138546 >>> max(0,strike-ST()) 12.364539627699976 >>> max(0,strike-ST()) 6.9642390525076365 >>> max(0,strike-ST()) 0 >>> max(0,strike-ST()) 8.417050607508223 >>> max(0,strike-ST()) 0 >>>
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 4.147540635574047 >>>
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 rate=.0535 vol=.34 delta=lambda x:x*exp((rate-.5*vol**2)*(1.0/365)+vol*(1.0/365)**.5*gauss(0,1))-x St0=41.75 time=61 ST=lambda: reduce(lambda x,y: x+delta(x),range(time),St0) strike=45 trials=10000 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.