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