Monte Carlo Application: Portfolio Diversification via Strategy Diversification

Tastytrade did a segment in May 2016 on strategy diversification, looking at the correlation between different option strategies and holding SPY stock outright.  This was the first time I thought of using options as a means of gaining some portfolio diversification.  In retrospect, this ability to diversify using options should have been obvious.  After all, a short premium trade can be profitable if the underlying doesn't move at all, or we can create a position that moves inversely to the direction of the underlying, or even have a positive correlation given a certain percentage move, but a negative given a larger.  But as I said, I had never considered this as a means of diversification.

Tastytrade used several options positions with SPY as the underlying.  We can use Monte Carlo methods to get at similar information.  We will use the assumptions used in the previous Monte Carlo articles.  Namely,

• Initial stock price = $100
• Annualized risk-free rate of return = 1%
• Implied volatility = 15%
• 30 trading days until option expiration
• 252 trading days per year

We will use the 95, 97.5, 100, 102.5 and 105 strikes in this example.  Just for the sake of naming conventions, we will round the option's delta to the nearest convenient number when tabulating and plotting the results.  Likewise, we will use the term "1 s.d. strangle" even though it's not quite one standard deviation.  Of course, the proper numbers are used internally.  We will use the following strategies listed in the table below.

Position Name	Strike(s)	Actual Delta
Long Stock	Not Applicable	1
Short ATM Call	100	52
Short 30 Delta Call	102.5	33
Short 20 Delta Call	105	19
Short ATM Put	100	48
Short ATM Put	97.5	-30
Short ATM Put	95	-15
Short Stradle	100	4
Narrow Strangle	97.5/102.5	3
1 s.d. Strangle	95/105	4
Long Call Spread	97.5/102.5	37
Long Put Spread	97.5/102.5	-37
Short Call Spread	102.5/105	-14
Short Put Spread	95/97.5	15

The stock Monte Carlo code is in the utilities.py file which I've included at the end of this post in the appendix.  After importing the needed libraries we write a simple function to calculate the option delta.  This is how the numbers in the table above were generated.  We do not need this code for the actual correlation computation but I left it in regardless.  We then define the parameters that go into the Monte Carlo calculation just as in previous articles.  We will also define all the strikes we will use.  Lastly, we run the Monte Carlo code.  It is more convenient this time not to have the results reshaped into a 2-dimensional matrix.  We will keep it as a 1-D vector.

from utilities import *
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats
import pandas as pd
import seaborn as sns

def delta(contract_type, d1):
    if contract_type.lower() == 'c':
        return scipy.stats.norm.cdf(d1)
    if contract_type.lower() == 'p':
        return scipy.stats.norm.cdf(d1) - 1
    
N_days = 30                     #  Number of trading days
N_sims = 1000                   #  Number of simulations

dpy = 252.0   #  Trading days per yeat
S0 = 100.   #  Initial stock price
K_atm = 100.0   #  At-the-money strike
K_call_105 = 105.0  #  Other strikes used
K_call_1025 = 102.5
K_put_975  = 97.5
K_put_95  = 95.0
sigma = 0.15   #  Anualized implied volatility
r = 0.01   #  Anuallized risk-free interest rate

#  Define the time until expiration.
t_days = np.arange(N_days, 0, -1)
t_days = np.append(t_days, 0.01)
t = t_days / dpy
t = np.tile(t, N_sims)  #  Duplicate and stack the time vector
    #  N_sims times

#  Our time step is one day
dt = 1.

#  Generate the Monte Carlo stock data
S = stock_monte_carlo(S0, N_days, N_sims, r, sigma, reshape = False)

Now that we have the stock prices, we'll calculate the values of $d_1, d_2$and the values of the calls and puts at the strikes we're using.

#  Functions called in this section imported from utilities.py

#  ATM calls and put

d1_atm, d2_atm = d(S, K_atm, r, sigma, t)
C_atm = call_price(d1_atm, d2_atm, S, K_atm, r, t)
d1_c_atm, d2_c_atm = d(S, K_atm, r, sigma, t)
P_atm = put_price(d1_atm, d2_atm, S, K_atm, r, t)
d1_p_atm, d2_p_atm = d(S, K_atm, r, sigma, t)

#  105 call

d1_c_105, d2_c_105 = d(S, K_call_105, r, sigma, t)
C_105 = call_price(d1_c_105, d2_c_105, S, K_call_105, r, t)

#  102.5 call and put
d1_c_1025, d2_c_1025 = d(S, K_call_1025, r, sigma, t)
C_1025 = call_price(d1_c_1025, d2_c_1025, S, K_call_1025, r, t)
P_1025 = put_price(d1_c_1025, d2_c_1025, S, K_call_1025, r, t)

#  97.5 call and put

d1_p_975, d2_p_975 = d(S, K_put_975, r, sigma, t)
P_975 = put_price(d1_p_975, d2_p_975, S, K_put_975, r, t)
C_975 = call_price(d1_p_975, d2_p_975, S, K_put_975, r, t)

#  95 put
d1_p_95, d2_p_95 = d(S, K_put_95, r, sigma, t)
P_95 = put_price(d1_p_95, d2_p_95, S, K_put_95, r, t)

Now that we have the individual option prices, we can calculate the prices of the spreads. When calculating the prices of strangles and straddles, I assumed long positions for the options and multiplied that by $-1$ for the actual correlation calculation.    For the vertical spreads, I worked out those details directly into the price of the spreads (ie, credit spreads have a negative price, debit spreads a positiveprice).  This should be obvious in the way the code is written.

#  Calculation of option spread prices

straddle = C_atm + P_atm
narrow_strangle = C_1025 + P_975
one_sd_strangle = C_105 + P_95
long_call_spread = C_975 - C_1025
long_put_spread = P_1025 - P_975
short_call_spread = C_105 - C_1025
short_put_spread = -P_975 + P_95

We're going to put the data in a Pandas data frame for analysis.  To that end, we'll first build a dictionary to hold the data.  The dictionary keys will end up being the labels in our frame and on the plot.

d =    {'Long stock': S, \
        'Short ATM call': -C_atm, \
        'Short 30 Delta call': -C_1025, \
        'Short 20 delta call': -C_105, \
        'Short ATM put': -P_atm, \
        'Short 30 Delta put': -P_975, \
        'Short 15 Delta put': -P_95,
        'Short straddle': -straddle, \
        'Narrow strangle': -narrow_strangle, \
        'One sd strangle': -one_sd_strangle, \
        'Long call spread': long_call_spread, \
        'Long put spread': long_put_spread, \
        'Short Call Spread': short_call_spread, \
        'Short put spread': short_put_spread}

#  Create the dataframe
df = pd.DataFrame(d)

#  Rearrange the columns

#  Place frame columns in the order I want
df = df[['Long stock',
        'Short ATM call',
        'Short 30 Delta call',
        'Short 20 delta call',
        'Short ATM put',
        'Short 30 Delta put',
        'Short 15 Delta put',
        'Short straddle',
        'Narrow strangle',
        'One sd strangle',
        'Long call spread',
        'Long put spread',
        'Short Call Spread',
        'Short put spread']]

Now we calculate the correlations and plot the results using a heatmap.

#  Calculate correlation dataframe
corr = df.corr()

#  Code for plotting the heatmap.  We include a mask since the matrix is
#  symmetrical.  This will get rid of the upper triangular portion.  We'll
#  include the diagonal though.
sns.set(style = 'white')
mask = np.zeros_like(corr, dtype = np.bool)
mask[np.triu_indices_from(mask)] = True
np.fill_diagonal(mask, False)

#  Create the figure window, set the color pallette, and create the plot
fig, ax = plt.subplots( figsize = (15,15) )

cmap = sns.diverging_palette(220, 10, as_cmap = True)

sns.heatmap(corr, mask = mask, cmap = cmap, annot = True, vmax = 1, center = 0,
        vmin = -1, square = True, linewidth = 0.5, cbar_kws = {'shrink': 0.5})

plt.show()

The resulting plot is shown below if Fig 1.The results agree with out intuition, a long call spread is strongly correlated with the stock movement, a strangle is relatively uncorrelated, etc.  They are also similar to those found by tastytrade.

Figure 1:  Heat map and correlation values among the various positions.

This calculation was straightforward and I think this will be the last post in this series on simple Monte Carlo techniques.   If I revisit Monte Carlo methods it will likely be applied to fields other than finance.

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 III:   Monte Carlo Techniques: Calculating the Probability of Making 50% of Max Profit on Short Option Positions

Appendix

The code for the Monte Carlo and options pricing model are contained in a separate file, utilities.py, shown below.

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

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