I was looking for a mathematical formula for the probability of a stock price touching a strike price of an option any time during it’s lifetime, assuming the stock price follows a traditional Black-Scholes model of Brownian Motion. I found some traction from from http://quant.stackexchange.com/questions/235/probability-of-touching, but unfortunately the equation did not seem to work according to freely available online calculators. My understanding is that ThinkOrSwim has this as a “proprietary” calculation on their options, however I don’t have that platform, and would like to calculate this “in house,” so to speak.

Anyway, one of the answers to the previous post linked to https://github.com/barrycarter/bcapps/blob/master/box-option-value.m, which is a Mathematica script for calculating the probability of a stock price hitting a “box” in a matter of x hours.

It follows the equation ..

\(\frac{{\partial}}{{\partial m }}P(m,t_1 ,t_2 ) =\frac{{1}}{{\sqrt {2 \pi ^2 t_1 (t_2 - t_1 ) } }}\int_{ - \infty }^m {\exp \bigg[ - \frac{{(m - x)^2 }}{{4(t_2 - t_1 )}} - \frac{{x^2 }}{{2t_1 }}\bigg]{\rm d}x}.\)

Now, I am not a big math person, but this is saying this is the probability of a particle in Brownian Motion hitting “m” between time “t1” and “t2”. The math here isn’t particularly important (well it actually is the most important, but for people like me who have to stand on the shoulders of someone who understands how to actually come to the above equation, it is just a piece of information in a puzzle).

So, after experimenting with the best way to actually use this script, I first converted it to Sage, an open source mathematics program that only works on Linux. This was a kind of cludgy solution that I did not like since my “toolbox” of finance programs that I build is run on .NET on Windows. Eventually I found that I could do the entire calculation in Python using the open source SciPy and NumPy. After a failed attempt at getting those to run on IronPython for .NET (which honestly should work, it is just in its infancy really and I think it is simply a bug), I decided to just use native Python on Windows. This works well for my purposes since the program can simply call the python script at the command line, pass the appropriate parameters, and then parse the result.

Here is the above script that was written in Mathematica, rewritten into Python using SciPy and NumPy.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

import sys
import scipy.stats
import scipy.integrate
import scipy.special
import math
from scipy import inf

if (len(sys.argv) < 7):
print "Must specify all variables"

p0 = float(sys.argv[1])
v = float(sys.argv[2])
p1 = float(sys.argv[3])
p2 = float(sys.argv[4])
t1 = float(sys.argv[5])
t2 = float(sys.argv[6])

def pdflp(x,p0,v,t1):
return scipy.stats.norm(math.log(p0), math.sqrt(t1)*v/math.sqrt(365.2425*24)).pdf(x)

def cdfmaxlp(x,v,t1,t2):
return (1-scipy.special.erf(x/(v*math.sqrt(t2-t1)/math.sqrt(24*365.2425/math.sqrt(2.0)))))

def upandin(p0,v,p1,p2,t1,t2):
integrand = lambda x: (pdflp(x,p0,v,t1) * cdfmaxlp(math.log(p1)-x,v,t1,t2))
return scipy.integrate.quad(integrand,-inf,math.log(p1))

def hitleftedge(p0,v,p1,p2,t1,t2):
integrand = lambda x: (pdflp(x,p0,v,t1)*1)
return scipy.integrate.quad(integrand, math.log(p1), math.log(p2))

def downandin(p0,v,p1,p2,t1,t2):
integrand = lambda x: (pdflp(x,p0,v,t1) * cdfmaxlp(x-math.log(p2),v,t1,t2))
return scipy.integrate.quad(integrand, math.log(p2), inf)

def probabilitytouch(p0,v,p1,p2,t1,t2):
return (upandin(p0,v,p1,p2,t1,t2) + hitleftedge(p0,v,p1,p2,t1,t2) + downandin(p0,v,p1,p2,t1,t2))

print probabilitytouch(p0,v,p1,p2,t1,t2)</pre>

Now you can just run this script from the commandline, like “python scriptname currentprice volatility strike1 strike2 time1 time2”. So, for example, let’s take SPY, the popular ETF for the S&P 500. Right now the price of the stock is $128.60, the volatility is 21.49%, and let’s look at the probability of hitting a strike of $130.00. Notice that we are not using multiple strikes to price a box like the original script intends, which is fine. And an expiration of November 19th, which is 21 days away (21*24 is 504 hours).

I would simply run the script in a command prompt (assuming python, SciPy, and NumPy are installed).

c:/python27/python.exe "c:/myscript.py" 128.60 .2149 130.0 130.0 1 504

This gives the output ..

(0.8027563748048097, 8.583401667801735e-09, 0.0, 0.0,
1.1850503070997317e-06, 1.1480513621859987e-11)

There are 6 numbers here separated by commas. All together they represent the probability of the price hitting the “box” between strike1 and strik2 (which in our case is simply a line because they are the same value, $130). Each probability has a number after it (such as 8.583401667801735e-09) which represents the error of the calculation, so values 1, 3, and 5 are the ones we are looking at. We add the first value, 0.8027563748048097, to the third value, 0.0, to the 5th value, 1.1850503070997317e-06. Now since the 5th value is so close to 0 we count it as 0. The probability is therefore .8027 or 80.27% of the stock hitting $130 between now and 21 days from now.

Keep in mind this follows the idea of Brownian Motion of stocks and all the rules of the Black-Scholes formula and is not “law” for where the stock will be 21 days from now. That is to say, this is a mathematical approximation of an approximation of a process that cannot really be described mathematically. So why is it valuable?

Well, the way I am using it is not to predict where a price will be, but where it has a low percentage to be, or more simply to find where the price has a high probability of not being on a certain date. This helps in a strategy where you are selling options (more specifically, veritcal credit spreads). We want to know the probability of a short strike being hit, and typically we want a 20% or lower chance of a short strike being hit. This formulation gives us something to quantitatively calculate and shoot for. It is important that this not take on the only piece of the decision making process though, and you should use some kind of outlook or technical analysis to further analyze where you expect the stock to go.