Binary Options and Implied Distributions with Python

by John | December 28, 2020



A binary option is a type of derivative in which a fixed payoff is received should the asset reach a certain level at expiration. A binary option with a payoff of 1 is often known as a digital option. These options are very similar to bets due to their relative simplicity. We can get some nice mathematics intuition regarding option pricing through studying binaries, which I hope to share with you today. 




In this article we will give an explanation of the mathematics behind binary option pricing along with a Python implementation for closed form and Monte Carlo pricing techniques. 

To finish off this article we will then give an example of getting the implied distribution of the stock price at expiration using binary options. 



It is worth mentioning at this point, that Binary options have been the subject of much controversy with regulators having worries about protecting investors from what is often outright fraud. Countries such as Canada, Germany and Israel have went as far as outright banning the sale of binary options to retail clients. In the United Kingdom, at one stage binary options were regulated by the Gambling Commission (FCA regulated now) hopefully this illustrates the point that the author does not recommend trading binary options unless serious due diligence is done. This article should be viewed as an educational resource as opposed to a promotion of trading these instruments for real money. A possible rule of thumb for discriminating between options providers is : Do they offer products that with an expiry of less than 1 minute? If yes, then it might be better to find another broker. 


 With that said let's begin! 


Simulation Method 


Consider an option that pays a fixed amount x conditional upon some event occurring in the market. Take an example of a stock currently trading at $100 with a binary option that pays $5 in the event the stock is greater than $115 in 3 month's time. Note that it doesn't matter whether the stock is $200 or $116 for an option of this nature, the payoff is $5 regardless.

You estimate the volatility to be 30% per annum. And the appropriate interest rate is 2% per annum. 


So the question is now how to price such as instrument?  

The reader may realize that it is useful to consider the question above as a probability question, in that we are asking how often would the stock finish above the strike. 


First we will calculate this by simulation as this is perhaps the most intuitive way to look at a problem of this nature. Below are the steps to complete this pricing method. Note we are assuming a log-normal distribution of stock prices at expiry, which is rather unrealistic but should serve to illustrates the concept. 


1) Simulate \(S_T\) according to Geometric Brownian Motion


See this article on where it comes from. Let N in the second line below be the number of draws to take from the distribution. 


\(ln(S_T) - ln(S_0) \sim \mathcal{N} (\ (\mu - \frac{\sigma^{2}}{2} ) T, \ \sigma^{2} \sqrt{T}) \)


\(S_T = S_0e^{(r- \frac{\sigma^2}{2} )T + \sigma dW_T^i} , \ for \ i \in [1,2 , ... ,N]\)



Below we simulate 10 million terminal stock prices, this should be sufficient to get a good approximation of the true distribution of stock prices at expiration. 


S=100.0 # spot stock price
K=115.0 # strike
T=0.25 # maturity 
r=0.02  # risk free rate 
sigma=0.3  # annualized volatility
Ndraws = 10_000_000

dS = np.random.normal((r - sigma**2/2)*T , sigma*np.sqrt(T), size=Ndraws)
ST = S0 * np.exp(dS)

n, bins, patches = plt.hist(ST,bins=250);
plt.title('Stock Simulation')



binary option simulation for Stock



 2) Calculate how often The stock is greater than the strike price. 

Imagine zipping along the x axis of the histogram above, and adding one to the total if the stock price from the draw is greater than the strike. We then count the number of ones and divide this sum by the number of draws which is 10 million in this case. 

The formula for this is shown below where \(1\) is the indicator function returning 1 or 0 dependent upon the condition being satisfied. The formula below represents the probability the stock is above the strike at expiration. Arguably we should we using an integral here as in the previous simulation but hopefully this way is more intuitive. 


\(\mathbb{P}(S_T >K) =\frac{1}{N} \displaystyle \sum_{i=1}^N \mathbb{1}_{\{ S_T > K\}} \)


The script below shows that the simulation approximates this probability as 16.5%. This should not be confused with the risk-neutral probability. Although viewing the formula here should give a good intuition as to what exactly a risk-neutral probability actually is when we encounter it later on in the article. 


n, bins, patches = plt.hist(ST,bins=250);
plt.title('Stock Simulation')
for c, p in zip(bins, patches):
    if c > K:
        plt.setp(p, 'facecolor', 'green')
        plt.setp(p, 'facecolor', 'blue')

timesInMoney = len(ST[ST>K])

P = timesInMoney/Ndraws
print('The Option was in the money ', f'{timesInMoney} times')
print(f'Out of {Ndraws} simulations')
print(f'Probability of being in the money {P}%')

#The Option was in the money  1651856 times
#Out of 10000000 simulations
#Probability of being in the money 0.1650933%


From the script above we see that the stock will be greater than the strike approximately 16.5% of the time. 


probability stock is above certain level at expiration



3) Calculate Value of Option.

Let \(Q\) be the payoff from the binary option ($5 in our example) The fair price for the binary option is then as follows.


\(\text{Binary Call Value} = \mathbb{P}(S_T>K) \times Qe^{-rT}\\ \\ 0.82=0.165\times 5e^{(-0.02 \times \frac{1}{4})} \)


So the fair price here is approximately $0.82. You can get the risk-neutral probability by multiplying \(\mathbb{P}(S_T>K) \times e^{-rT}\) if you so wish.


The same procedure can be applied to put options by changing  \(\mathbb{P}(S_T>K) \times e^{-rT} \ \text{to} \ \mathbb{P}(S_T<K) \times e^{-rT}\) the function below can be used to try different options parameters. 


def monte_carlo_binary(S, K, T, r, sigma, Q, 
                       type_='call', Ndraws=10_000_000, seed=0):
    dS = np.random.normal((r-sigma**2/2)*T, sigma*np.sqrt(T),size=Ndraws)
    ST = S *np.exp(dS) 
    if type_ =='call':
        return len(ST[ST>K])/Ndraws * Q *np.exp(-r*T)
    elif type_ == 'put':
        return len(ST[ST<K])/Ndraws *Q*np.exp(-r*T)
        raise ValueError('Type must be put or call')

monte_carlo_binary(S, K, T, r, sigma, 5)  





Black-Scholes Closed Form

We can also use the Black-Scholes formula to price binary options, for this we will need the d2 from the previous article.  The formulae for calls and puts are given below. 


Call formula and Python Implementation




from scipy.stats import norm

S=100.0 # spot stock price
K=115.0 # strike
T=0.25 # maturity 
r=0.02  # risk free rate 
sigma=0.3  # annualized volatility

def d2(S, K, T, r, sigma):
    return (np.log(S/K) + (r- sigma**2 / 2)*T) /  (sigma*np.sqrt(T))

def binary_call(S, K, T, r, sigma, Q=1):
    N = norm.cdf
    return np.exp(-r*T)* N(d2(S,K,T,r,sigma)) *Q


Let's just take a moment to equate some concepts from the Monte-Carlo method we discussed. Notice that when we pass Q=1 in we get a closed form solution to the probability the stock will finish in the money. 


\(\mathbb{P}(S_T>K) \times e^{-rT} = N(d1)\)


\(0.165\times e^{(-0.02 \times \frac{1}{4})} = N(d1)\)


binary_call(S, K, T,  r , sigma, Q = 1)



Notice that we can recover the probability value we got from the Monte-Carlo simulation by the following. 


\(\mathbb{P}(S_T>K) = N(d1)\times e^{rT} \)


\(0.165 = e^{(0.02 \times \frac{1}{4})} N(d1)\)


binary_call(S, K, T,  r , sigma, Q = 1) * np.exp(r*T)



And Pricing our example option we get approximately the same value. Increasing the Ndraws parameter will reduce this error, however we see below it is fairly accurate and they are in fact measuring the same quantity. 


binary_call(S, K, T,  r , sigma, Q =5)

monte_carlo_binary(S, K, T, r, sigma, 5)  



Put formula and Python Implementation


The formula for pricing a binary put option is given below, in this case we are measuring the probability of the stock being below the strike price. 




def binary_put(S, K, T, r, sigma, Q=1):
    N = norm.cdf
    return np.exp(-r*T)* N(-d2(S,K,T,r,sigma)) *Q


Let's try that formula out on pricing a put option with the same parameters as the call we have used throughout this article


binary_put(S, K, T,  r , sigma, Q = 1)

monte_carlo_binary(S, K, T, r, sigma, 5, 'put')
# 4.153253729048758


Now consider if we could have inferred this value without actually using either formula. Since we know that the problem is binary i.e. one of the two events must occur, the stock is either above the strike or below it, the following relationship must hold


\(\mathbb{P}(S_T>K) + \mathbb{P}(S_T<K) = 1\)


To adjust this for a risk neutrality argument we can state the equality shown below. Let \(C_k\ \text{and} \ P_k\) be the prices of digital calls and puts respectively. 


\(e^{rT}(C_k + P_k) = 1\)


np.exp(r*T)*(binary_call(S,K,T, r, sigma, Q=1)     +binary_put(S,K,T, r, sigma, Q=1))



And then it follows that to price a non-digital binary:


\(e^{rT}(C_K + P_K) = Q\)



Clearly once we know the price of a binary call option we can then infer the price of the put. 


\(P_K = e^{-rT}Q -C_K\)


Q = 5
print(np.exp(-r*T)*Q - (binary_call(S,K,T, r, sigma, Q=Q)) )
# 4.153311176925953

print(binary_put(S,K,T, r, sigma, Q=Q))






The following points are important to remember in order to fully grasp the next section. 


1) \(N(d1)\)  from the Black-Scholes formula represents the risk neutral probability of the stock being above the strike at expiration. 

2) \(N(-d_2)\) is the risk neutral probability of the stock finishing below the strike at expiration. 





Implied Probability Distribution from Market Data


In this mini project we will take some of the things we have learned about binary options and apply them to some real market data. It may be useful to read this article on implied volatility if you are unfamiliar with the concept. 


The goal of this section is to create a cdf and pdf of the market's expectations regarding the price of Apple stock on the 19th of February.


To follow along you can either download the market data yourself from github here or you can simply download it using Pandas as shown below. 


import pandas as pd

df = pd.read_csv('')


#    K        IV       S         T
#0  55  1.017583  131.97  0.207843
#1  60  0.923829  131.97  0.207843
#2  65  0.871095  131.97  0.207843
#3  70  0.802736  131.97  0.207843
#4  75  0.736331  131.97  0.207843


Notes on the data

1) The data was collected from Yahoo Finance prior to the market open on 28th December 2020 


2) All of the options expire on the 19th February 2021


3) The Time to expiration variable T was calculated by taking the number of days and dividing by 255.(Could be more accurate admittedly) 


4) It is not clear which value Yahoo Finance uses to calculate implied volatility, however, we won't be dealing with market prices and therefore are making some unrealistic assumptions in order to illustrate the concept. Feel free to try it on different data


5) We will assume the risk free rate is 0.1% per annum taken from here






Interpolate and Extrapolate Implied Volatility 


Here we use a polynomial fit with degree 5 to get our new implied volatility values. Since the highest and lowest strike available is 230 and 55 respectively we are going to extrapolate for values between 1 - 55 . While we do suspect that values towards the end of this distribution are highly likely to be much higher in real life, we will use the following model simply for illustrative purposes. 


vols = df.IV.values #volatility values
Ks = df.K.values #Strikes 

newK = np.arange(1, 270,0.0001) # new higher resolution strikes

poly = np.polyfit(Ks,vols,5) #create implied vol fit
newVols = np.poly1d(poly)(newK) # fit model to new higher resolution strikes

plt.plot(newK, newVols)
plt.title('Implied Volatility Function')
plt.ylabel('Implied Vol')


implied vol interpolation


So what we have now is a method to approximate the appropriate volatility values from the data we collected from Yahoo Finance. The reader is encouraged to play around with the function below and compare it with the plot above. 


def vol_by_strike(polymdl, K):
    return np.poly1d(polymdl)(K)

vol_by_strike(poly, 131.97)


So we see from above that the dead on the money volatility is 40.39% 



Create Risk Neutral Cumulative Distribution Function for Stock Price at Expiration. 


To create a cdf we will want to calculate the weight to the left of the given point, the aforementioned point here is the strike. Referring back to the examples at the beginning of the document we know to calculate this value we can use a digital put option. 


S = df.S[0] # extract S_0
T = df['T'][0] #extract T
r= 0.001 # Chosen for convenience

binaries = binary_put(S, newK, T, r, newVols) #calculating the binaries

plt.plot(newK, binaries, label='cdf')
plt.axvline(S,color='black',linestyle='--', label='$S_0$')
plt.title('Cdf of Apple on 19th Feb 2021')


apple cdf implied by calls



We can then use the following script to assign a probability to each value of \(S_T\) , I hesitate to call this a pdf since we are assigning discrete probabilities to individual price. However, it is useful for illustrative purposes. We will also add a constant volatility distribution i.e. the implied vol at \(S_0\) for comparative purposes. 


binaries = binary_put(S, newK, T, r, newVols)

PsT = []  # market implied probabilities

for i in range(1, len(binaries)):
    p = binaries[i] -binaries[i-1]

constant_vol = vol_by_strike(poly, S)
binaries_const = binary_put(S, newK, T, r, constant_vol)

const_p = [] 

for i in range(1, len(binaries_const)):
    p = binaries_const[i] -binaries_const[i-1]

plt.plot(newK[1:], PsT, color='red',label='Market Implied')
plt.plot(newK[1:], const_p, color='black',
         label='Constant Vol')
plt.title('Market implied vs Constant')


Notice that using a constant volatility, there is effectively 0% probability that the stock is below 60 at expiration. However, the market doesn't agree with this idea, perhaps we can interpret this as the risk rare events such as war , natural disaster etc. 


implied pdf by calls vs constant vol model



Let's explore what we can do with this distribution now that we have it. Let's see how we can calculate the probability that the stock is within a certain interval on the expiration date. 

For example say we want to find that risk-neutral probability of :


\(\mathbb{P}(110 < S_T < 140) \)


PsT =np.array(PsT)
idx = np.argwhere((newK[1:] <= 140) & (newK[1:] > 110))

# 0.4409622114337626


So according to the market there is a 44.10% chance of the Apple expiring in this interval illustrated below.


probability of stock finishing within range of prices



plt.figure(figsize=(7, 4))
plt.plot(newK[1:], PsT, color='black')
plt.fill_between(newK[1:], PsT, where=((newK[1:] < 140)& (newK[1:] > 125)),
                 facecolor='black',alpha=0.9 )
plt.title('Probability  Interval')
plt.xticks(np.arange(0,230, 20.0));




Trading Strategy Implications

Recall the strategies illustrated in previous articles here and here. Let's say a trader has his own model for estimating the density of \(S_T\) at expiration, if the model predicts a larger density than the market is implying, the trader could take a long position in a butterfly spread or if the opposite were true he could take a short position in a butterfly. 



Hopefully this article has helped you make a connection between probabilities implied by option prices and also an intuitive understanding of risk-neutral probabilities and what they actually mean. 



Join the discussion

Share this post with your friends!