Monday, January 22, 2018

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

No comments:

Post a Comment