Monte Carlo Techniques: Calculating the Probability of Making 50% of Max Profit on Short Option Positions

Short options positions--  for the time being, we'll look only at naked positions--   have a defined, finite level of profitability and potentially unlimited loss.  Despite the large potential losses, these trades have a high probability of success.   The folks at tastytrade have shown managing these trades at 50% of the maximum profit gives better performance than simply holding to expiration.  Subsequent research by them over the years have corroborated their initial finding.

This probability of a short option expiring in-the-money can be calculated from the Black-Scholes option pricing model.   This is quite straightforward.  But there is no simple way to calculate the probability of making 50% of max profit on a trade or for that matter estimating the probability of at some point in the trade's lifetime, being down a certain multiple of the premium collected.

Tastyworks, a relatively new brokerage, displays the probability of making 50% on options trades where applicable.  They use a Monte Carlo method to arrive at this number.  In this article, we will use the ideas from the previous article on modeling stock price action to arrive at a similar number.  As before, we'll be doing this in Python 2.7.

Recall from my article on calculating implied volatility that the price of a call $C$ or put $P$ is given by the Black-Scholes formula,
$$C = \Phi(d_1) S - \Phi(d_2) K e^{-r t},$$ and
$$P =  \Phi(-d_2) K e^{-r t} - \Phi(-d_1) S, $$ respectively.

The terms $d_1$ and $d_2$ are given by,
$$d_1 = \frac{1}{\sigma \sqrt{t}} \left[ \ln\left(\frac{S}{K}\right) + \left(r + \frac{\sigma.^2}{2}\right) t\right],$$
and
$$d_2 = d_1 - \sigma \sqrt{t}.$$

In the above equations, $S$ is the stock price, $K$ is the strike price, $\sigma$ is the implied volatility,  $r$ is the annualized risk-free interest rate, $t$ is the time remaining until expiration denoted in years, and $\Phi$ is the normal cumulative distribution function.

For this example, we'll use a call option.  The same logic applies to puts.  We will use the following values in our calculation.

• $\sigma$ = 0.15
• $r$ = 0.01
• $t$ = 30 days
• Initial stock price, $S_0$ = $100
• Strike price $K$ = $105

Let's put the code for the stock Monte Carlo calculation into its own file along with functions to calculate $d_1$, $d_2$, and the prices of calls and puts.  Then we can import this file to use as a library in later code.  The only change I've made to the stock_monte_carlo function from the previous article is to have an option that turns off reshaping the results into a matrix where each row corresponds to a simulation and each column is that day's price.  In the next article, it will be  more convenient to have all the results returned as a single vector, so I've included that option in the code below.

#  utilities.py
import numpy as np
from math import sqrt
import scipy.stats
import scipy.sparse
import scipy.sparse.linalg

def stock_monte_carlo(init_price, N_days, N_sims, r, sigma, reshape = True):

    #  Scale interest rates and volatility.  Define time-step.
    r = r / 252.
    dt = 1.0
    sigma = sigma / sqrt(252.0)

    #  Calculate vector of normally distributed numbers and use it to
    #  calculate the daily percent change.
    epsilon = np.random.normal( size = (N_sims * N_days + N_sims - 1) )
    ds_s = r * dt + sigma * sqrt(dt) * epsilon

    #  Step up matrix diagonals
    ones = -np.ones( (N_sims * N_days + N_sims) )
    ones[0:-1:N_days+1] = 1.

    ds_s[N_days:N_days * N_sims + N_sims:N_days+1] = -1
    d = [ds_s + 1, ones]
    K = [-1, 0]

    #  Solve the system of equations
    M = scipy.sparse.diags(d, K, format = 'csc')
    p = np.zeros( (N_sims * N_days + N_sims, 1) )
    p[0:-1:N_days+1] = init_price
    s = scipy.sparse.linalg.spsolve(M, p)

#  Reshape the column vector so the function returns a matrix where
#  each row is a single simulation with each day corresponding the the
#  columns.  This is dine by default but can be overridden by the user
    if reshape == True:
        s =  np.reshape(s, (N_sims, N_days+1))

    return s

def call_price(d1, d2, S, K, r, t):
    C = np.multiply(S, scipy.stats.norm.cdf(d1)) - \
    np.multiply(scipy.stats.norm.cdf(d2) * K, np.exp(-r * t))
    return C

def put_price(d1, d2, S, K, r, t):
    P = -np.multiply(S, scipy.stats.norm.cdf(-d1)) + \
    np.multiply(scipy.stats.norm.cdf(-d2) * K, np.exp(-r * t))
    return P

def d(S, K, r, sigma, t):
    d1 = np.multiply( 1. / sigma * np.divide(1., np.sqrt(t)), 
        np.log(S/K) + (r + sigma**2 / 2.) * t  )
    d2 = d1 - sigma * np.sqrt(t)
    return d1, d2

As mentioned above, this is actually a very straightforward calculation.  As in the last article, we will start by importing the needed libraries and set the random seed to a constant for debugging purposes.  We will also turn off numpy's divide-by-zero warning.  This is not a good practice in production code, but will work for us here.  The issue is that there is a $1/t$ term in the calculation for $d_1$ and when there is zero time to expiration, that term goes to $-\infty$.  Numpy takes $x/0$ to be $\infty$ for nonzero values of $x$. Since $\Phi(-\infty)$ is zero, this will cause an out-of-the-money option to have zero value at expiration and an in-the-money option to have only intrinsic value.  Again, this works for us, but is not a good practice in general.  We will also set the random seed to a constant for debugging purposes.

from utilities import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats

np.random.seed(2)
np.seterr(divide = 'ignore')

Next, just as in the last post, we set up our parameters and run the code to generate the needed stock data.

N_days = 30                     #  Number of trading days
N_sims = 1000                   #  Number of simulations

#  Set up our parameters
dpy = 252.0                             #  Trading days per year
S0 = 100.                               #  Initial stock price
K = 105.0;                              #  Strike price
sigma = 0.15                            #  Implied volatility
r = 0.01                                #  Risk free rate

#  Since we need te price of the call option at every day, we'll generate a
#  sequence days from N_days to zero
t_days = np.arange(N_days, 0, -1)
t_days = np.append(t_days, 0.)
t = t_days / dpy                        #  Convert days to yeears

dt = 1.                                 #  Time step in days

#  Run the Monte Carlo code to generate stock data
S = stock_monte_carlo(S0, N_days, N_sims, r, sigma)

Now let's calculate the vales of $d_1$, $d_2$, and the call prices $C$.  Then we plot all the call prices.

#  Calculate the values of d1, d2, and the call Prices
d1, d2 = d(S, K, r, sigma, t)
C = call_price(d1, d2, S, K, r, t)

#  Plot calll prices
#t = np.arange(0, N_days + 1)
for i in range(N_sims):
    plt.plot(t * 252, C[i,:], 'k', alpha = 0.05)

plt.xlabel('Time to Expiration (days)')
plt.ylabel('Call Price ($)')
plt.grid(True)
plt.gca().invert_xaxis()
plt.show()

As with the last article, we will make the lines translucent so we can get a good visualization of the results.   The transparency gives us a decent indication of the probability of those prices being realized.  The darker the region, the more likely the result.

Figure 1:  Simulated call prices  plotted as a function of time until expiration.

As we can see, it looks like the majority of the calls eventually expire worthless.  Recall that we are short the call, so this is a good result for us.  There are however, some outlier moves to the upside which would result in large losses.  Let's visualize this using a histogram instead.

plt.figure()
plt.hist( C[:,-1], bins = 25, edgecolor = 'k', linewidth = 0.5 )
plt.xlabel('Call Price ($)')
plt.ylabel('Count' )
plt.yscale('log', nonposy='clip')
plt.show()

The histogram is plotted below in Fig 2.  Since almost all calls go out worthless or nearly so, we need to plot the Y-axis on a log scale.

Figure 2:  Histogram showing the distribution of call prices for 1000 simulated trades.  Since the overwhelming majority expire worthless or nearly so, we plot the counts on a log scale.

I find it a good practice to build "sanity checks" into these types of calculations. Let's do the following:

1. Calculate the total number of winning and losing trades
2. Make sure the number of winners plus the number of losers adds up to $N_\mbox{sims}$..
3. Use the information above to calculate the probability of profit at expiration based on the Monte Carlo calculations.
4. Compare this with the theoretical probability of profit from the Black-Scholes model.

#  What was the initial premium collected for the sale of the call?
initial_price = C[0,0]
print 'Initial call price = ', initial_price

#  P&L if held to expiration.  At expiration,  many losing trades are there?
losers = C[:, -1] > initial_price
print 'Number of losing trades:  ', np.sum(losers)

winners = C[:, -1] <= initial_price
print 'Number of winning trades:  ', np.sum(winners)

#  Sanity check by making sure the number of winners + losers equals the
#  total number of simulated trades
print 'Sanity check: ', np.sum(winners) + np.sum(losers), ' == ', N_sims
print '\n'


d1, d2 = d(100., K + initial_price, r, sigma, t[0])
print 'Probability of profit from Black-Scholes:  ', 1 - scipy.stats.norm.cdf(d2)

print 'Estimated probability of profit from Monte Carlo:  ', \
    np.sum(winners) / float(N_sims)

Running this produces the results shown in Fig. 3.

Figure 3:  Results of some simple calculations to see if numbers are reasonable.

We see these numbers make sense.  The number of losers plus the number of winners does add up to the total number of simulated trades, and the probability of profit from the Monte Carlo Simulation agrees nicely with the number obtained from Black-Scholes.

Let's make one more check to see if these numbers are reasonable before we move on. We are assuming that the implied volatility matches the actual realized movement of the stock, and that the volatility remains constant over the life of the trade. Both of these assumptions are false, but if we take them as true, we'd expect to be playing what amounts to a zero sum game. Over time, our profits and losses should cancel out. Let's see how accurate that is.

The following code snippet finds the indices of the winning trades, then calculates the some of profits from all winning trades.  Then we repeat the process for the losers.

win_ind = np.where( C[:, -1] <= initial_price )
wins = C[win_ind, -1]
print 'Total profit from winners: ', np.sum( initial_price - wins )

loss_ind = np.where( C[:, -1] > initial_price )
losses = C[loss_ind, -1]
print 'Total losses: ', np.sum(losses - initial_price)

Running the above gives the results in Fig. 4.

Figure 4:  results of the code totaling all wins and losses.

As we can see, the numbers are pretty much equal. In this case, the win total is slightly more than the losses, but changing the random seed can change that result.  In any cases, the two numbers are always close to each other as wed'd expect.

OK, now that everything looks good, let's calculate the probability of making 50% of maximum profit on this trade.  To do this, will first look at every simulated price and see if it is less that or equal to half the collected premium.  The results of this operation will be a boolean matrix.

When the function np.sum is run on a boolean array, it interprets True as1 and False as zero.  We will sum across each trade.  If the trade is never at or above 50%, the result will be zero.  The result of this operation should be a single vector.  Then we just count the number of nonzero entries in this vector and divide it by the total number of trades.

half_max = initial_price / 2.
reached_half_max = C <= half_max
reached_half_max = np.sum(reached_half_max, axis = 1)
print 'Percentage of trades that reached 50% of max profit:  ', \
    np.sum( reached_half_max > 0 ) / float(N_sims)

This gives the following result.

There we have it  the percentage of making 50% of max profit on this this trades is roughly 92%.

We can turn this problem around.  90+% probability of profit on this looks really good.  Let's turn it around and see how much pain we are likely to endure in doing this.  What is the probability at some point, we will be a loser, down 100% the amount of premium collected.  In this example, we collect about $\$0.51$, so how likely is it that the call marks at $\$1.02$? We do it in the same way as above.

twice_max = initial_price * 2.
reached_twice_max = C >= twice_max
reached_twice_max = np.sum(reached_twice_max, axis = 1)
print 'Probability of loss at some point is 100% of the premium collected:  ', \
    np.sum( reached_twice_max > 0 ) / float(N_sims)

This gives the following result.

So while the probability of making 50% of max is excellent, there is about a 40% chance you'll be down some money on the trade at some point.  That isn't terribly surprising as there is a good chance our strike will be touched at some time throughout the trade.

In this article, we built upon the Monte Carlo code from the previous post in this series that simulates the statistical behavior of a stock.   We did the heavy mathematical lifting in that last post, and  bulk of this article boils down to taking those results and plugging them into the Black-Scholes model and analyzing the results.  The general concept here is straightforward and very powerful.  Though we've talked only about single short options, the same concept could be applied to long options and spreads of various types.  In the next article, we will look at the relations between options and their underlying stocks in terms of correlations.  Portfolio diversification is an important topic, and we can use options, perhaps somewhat counter-intuitively,  to achieve some of that diversification.

Articles in this series

Part I:  Monte Carlo Techniques: Estimating the Value of Pi with Random Numbers
Part II: Monte Carlo Techniques: Modeling Stock Price Action
Part IV:  Monte Carlo Techniques:  Estimating Correlations between Option Positions and  and Their Underlying Stock

Calculating Implied Volatility from an Option Price

Nowadays, with the ubiquity of computers and information, many traders wish to know how the sausage is made and work out for themselves some of the numbers appearing on the screen of their trading platforms.  This is often a simple matter of looking up the formulas involved and plugging in the numbers.

When calculating the numbers pertaining to options, however, we run into an issue with implied volatility.  The Black-Scholes model tells us what an option should be worth given its strike price, the risk-free interest rate, the remaining time until expiration, the stock's price, and the implied volatility.  For example, the price of a call, $C$, is given by,

$$C = \Phi(d_1) S - \Phi(d_2) K e^{-r t},$$

where $S$ is the price of the stock, $K$ is the strike price, $r$ is the annualized risk-free rate, and $t$ is the remaining time to expiration expressed in years.  $\Phi$ is the normal cumulative distribution function, and $d_1$ and $d_2$ are given by,
$$d_1 = \frac{1}{\sigma \sqrt{t}} \left[ \ln\left(\frac{S}{K}\right) + \left(r + \frac{\sigma.^2}{2}\right) t\right],$$
and
$$d_2 = d_1 - \sigma \sqrt{t},$$
respectively.  In the above expressions for $d_1$ and $d_2$,  $\sigma$ is the implied volatility.

If we wish to calculate $\sigma$, we run into an issue.  All the variables in the above equations are known, and we can get the call price directly from the option chain, but we are incapable of isolating $\sigma$ algebraically.  Instead, we will have to turn to numerical methods to calculate the implied volatility.

Recall from high school algebra that if $y = f(x)$, the value of $x$ for which $f(x) = 0$ is called the root of the function $f$.  We will make use of a root finding algorithm to find our volatility, $\sigma$.  Our function will be the theoretical call price from the Black-Scholes model minus the known option price.  We will insist that equals zero and find the value of $\sigma$ that makes it so.  Explicitly, we want,

$$C(\sigma) - C_0 = 0,$$

where $C_0$ is the call price from the option chain.

For the sake if this article, we will assume the underlying stock is trading for 100, the strike price is 105, there are 30 days until the contract expires, and the risk-free rate is 1%.  We will plot $C(\sigma) - C_0$ for values of implied volatility from zero to one below just to get a sense of how this function behaves.

Figure 1.  We wish to find the value of $\sigma$ where $C(\sigma) - C_0$ is zero.

We wish to find the value of $\sigma$ where the function $C(\sigma) - C_0 = 0$.  From Fig 1., we can see this happens around $\sigma = 0.38$.  We will use an iterative method for finding an approximation to the root developed by Isaac Newton and Joseph Raphson.  The method requires an initial guess to start the process.  We will denote our successive approximations of the root as $\sigma_i, i = 0,1,2,...$ and will start with a guess of $\sigma_0 = 0.5$.

The Newton-Raphson method calculates the slope of the tangent line evaluated at the first iteration point.  It then uses the point where that tangent line crosses zero as the starting point for the next iteration.  This is shown graphically in Fig. 2.

Figure 2.  Plot of the line tangent to our function at $\sigma = 0.5$.  The intercept of this line with the x-axis is used as the starting point for the next iteration.  We can see after only one iteration it is relatively close to the actual root

The process then repeats until our test value for $\sigma$ is as close to zero as we like.  Fig. 3 shows the second iteration.  Note that after only two iterations, we are very close to the root.

Figure 3.  Tangent line for the second iteration

To implement this technique, we need to be able to calculate the derivative of the function of interest.  One of the nice features of the Black-Scholes model is that it has a closed form that can be easily differentiated.  Traders refer to the derivative of the option price with respect to volatility as vega and denote it with the Greek letter $\nu$.  Vega is given by,
$$\nu = \frac{\partial C}{\partial\sigma} = S \phi(d_1)\sqrt{t},$$
where $\phi$ is the normal probability density function.  With the above equations, we have enough information to implement a program to calculate the implied volatility of an option.  We will use Python for this exercise because it is a popular, freely available programming language that has a fairly extensive math and statistics libraries.

We will make use of the scypi.stats library as well as specific functions in the math package.  After importing these, we need to write a couple of helper functions to implement the Black-Scholes model for call option prices.  The first function calculates the values of $d_1$ and $d_2$.  The second function actually calculates the theoretical price of the call.

from scipy.stats import norm
from math import sqrt, exp, log, pi

def d(sigma, S, K, r, t):
    d1 = 1 / (sigma * sqrt(t)) * ( log(S/K) + (r + sigma**2/2) * t)
    d2 = d1 - sigma * sqrt(t)
    return d1, d2

def call_price(sigma, S, K, r, t, d1, d2):
    C = norm.cdf(d1) * S - norm.cdf(d2) * K * exp(-r * t)
    return C

Next we enter in the known option information as well as our initial guess for $\sigma$.

#  S  = Stock price
#  K  = strike
#  C  = price of call as predicted by Black-Scholes model
#  r  = risk-free interest rate
#  t  = time to expiration expressed in years
#  C0 = price of call option from option chain

S = 100.0
K = 105.0
r = 0.01
t = 30.0/365
C0 = 2.30

#  We need a starting guess for the implied volatility.  We chose 0.5
#  arbitrarily.
vol = 0.5

Next we define some variables needed for bookkeeping.  These will include a variable to count the number of iterations along with a maximum number of iterations.  We will use these to make sure our program doesn't get stuck in an infinite loop.  We will also set a tolerance where we will terminate the iterations when our test root gives a value that is close enough to zero to satisfy our needs.

epsilon = 1.0          #  Define variable to check stopping conditions
abstol = 1e-4          #  Stop calculation when abs(epsilon) < this number

i = 0                  #  Variable to count number of iterations
max_iter = 1e3         #  Max number of iterations before aborting

Lastly, we will do the iteration and print out the results.

while epsilon > abstol:
    #  if-statement to avoid getting stuck in an infinite loop.
    if i > max_iter:
        print 'Program failed to find a root.  Exiting.'

        break

    i = i + 1
    orig = vol
    d1, d2 = d(vol, S, K, r, t)
    function_value = call_price(vol, S, K, r, t, d1, d2) - C0
    vega = S * norm.pdf(d1) * sqrt(t)
    vol = -function_value/vega + vol
    epsilon = abs(function_value)

print 'Implied volatility = ',  vol
print 'Code required', i, 'iterations.'

Running this produces the result,

Implied volatility =  0.368856324914
Code required 3 iterations.

We see the program converges quickly taking only three iterations to find a value within our tolerance.  The above technique can be used for puts by substituting in the Black-Scholes formula for put prices, of course.   So that's it.  Pretty simple.

Newton's method works well when the function and its derivative are well-behaved.  In our case, the function only has one root, so we don't have to worry about the possibility of having several mathematically allowable answers and setting up our code and initial guess to converge to the proper root or deciding which of several results is the right answer to our problem.

Here is the above code in its entirety.

#!/usr/bin/python

from scipy.stats import norm
from math import sqrt, exp, log, pi

def d(sigma, S, K, r, t):
    d1 = 1 / (sigma * sqrt(t)) * ( log(S/K) + (r + sigma**2/2) * t)
    d2 = d1 - sigma * sqrt(t)
    return d1, d2

def call_price(sigma, S, K, r, t, d1, d2):
    C = norm.cdf(d1) * S - norm.cdf(d2) * K * exp(-r * t)
    return C

#  S  = spot
#  K  = strike
#  C  = price of call as predicted by Black-Scholes model
#  r  = risk-free interest rate
#  t  = time to expiration expressed in years
#  C0 = price of call option from option chain

S = 100.0
K = 105.0
r = 0.01
t = 30.0/365
C0 = 2.30

#  We need a starting guess for the implied volatility.  We chose 0.5
#  arbitrarily.
vol = 0.5

epsilon = 1.0  #  Define variable to check stopping conditions
abstol = 1e-4  #  Stop calculation when abs(epsilon) < this number

i = 0   #  Variable to count number of iterations
max_iter = 1e3  #  Max number of iterations before aborting

while epsilon > abstol:
    #  if-statement to avoid getting stuck in an infinite loop.
    if i > max_iter:
        break

    i = i + 1
    orig = vol
    d1, d2 = d(vol, S, K, r, t)
    function_value = call_price(vol, S, K, r, t, d1, d2) - C0
    vega = S * norm.pdf(d1) * sqrt(t)
    vol = -function_value/vega + vol
    epsilon = abs(function_value)

print 'Implied volatility = ',  vol
print 'Code required', i, 'iterations.'

For further reading see:

Wikipedia article on the Black-Scholes model
Wikipedia article on Newton's method

Updated 11/12/2018:  Fixed text which said the strike price was 155 with a stock price of 150.  The code uses a strike of 105 and a price of 100.  Sorry for not catching this sooner.