Options’ Implied Probability (with code): Risk-Neutral Densities modelling – Covid example
In this article, I will describe the method of retrieving the options’ implied probability from option prices. The result will be probability density functions, also referred to as “Risk-Neutral densities” (RND) or “State-Price densities” (SPD).
This article is a supplement to the previous one (check it out here). However, as many people were interested in getting the code, I am going to expand the topic and provide a code walkthrough. I won’t cover all the coding steps (running functions, etc.), but it is complete and ready to run.
Moreover, I am going to present you with a way to retrieve data from Theta Data, a very cool data provider, where you have institutional-quality data with retail prices. It will save us a ton of time, as getting rid of the bad and inaccurate data is a standard procedure when using common, open sources like Yahoo Finance.
Furthermore, as we have no limits to the data, I will show you RNDs from a very interesting time, which is the COVID breakout. It is approximately late February 2020. More about it later.
Why RND?
Options drive the market. They are heavily traded and a wonderful tool to hedge or trade (or go broke). These are probably something closest to the market expectations of the future today, and RNDs are an “implied view” of the future.
The RND is another tool that can help us comprehend the market and its prices. It can be useful to generate signals. I won’t delve into that today, but it is surely an interesting topic to explore.
RNDs act the same as regular density. There are multiple methods to retrieve them — parametric and nonparametric ones. Parametric methods depend on paramaters and have some predefined shape, that is influenced by parameters. Nonparametric means that there are no assumptions about the shape of the distribution or its parameters.
As it is difficult to have a garanular option chain for expirations, this entails that we will have to resort to some parametric solution. As I am skeptical about parameterizing RND, I will parametrize volatility instead. There are a ton of methodologies that are used in practice, such as the SABR model or others. We will use SABR, as there is a free library in Python, and then use the Breeden-Litzenberger formula. If something isn’t clear now, hold on tight.Data Acquisition
When delving into the topic of RNDs, the first crucial step lies in acquiring and refining the data. I will use Theta Data, as it is the best source out there. They have tons of data of the highest quality at affordable prices (check them out; the link is here). I highly recommend them.
I will use a standard monthly expiry in March (2020–03–20), so we can observe the market crash and the RND behavior. These are EOD quotes with the best asks and bids at that time. I am stressing that, very often, the EOD quotes can be last traded prices. It means that the price curve of the option can be jagged and unreliable. For example, the 400 strike could be traded at 15:55 and the 450 strike at 15:59, which will surely entail discrepancies around those prices. Nevertheless, I will show you a cool trick on how to handle this problem as well. The data is for the last 100 days before expiration. Here is the whole code showcasing how to fetch the data:
client = ThetaClient(username="yourmail@", passwd="your password")
def wrangle_response(response):
columns = response['header']['format'] + ['root', 'expiration', 'strike', 'right']
all_rows = []
for contract in response['response']:
contract_info = contract['contract']
for tick in contract['ticks']:
row = tick + [contract_info['root'], contract_info['expiration'], contract_info['strike'], contract_info['right']]
all_rows.append(row)
df = pd.DataFrame(all_rows, columns=columns)
return df
### If there are more expirations
def get_bulk_at_time_quotes(expirations, root):
url = "http://127.0.0.1:25510/bulk_at_time/option/quote"
headers = {"Accept": "application/json"}
all_data = pd.DataFrame()
successful_fetches = 0
for expiration in expirations:
print(f"Expiration download: {expiration}")
try:
start_date = expiration - timedelta(days=100)
formatted_expiration = expiration.strftime('%Y%m%d')
formatted_start_date = start_date.strftime('%Y%m%d')
querystring = {
"root": root,
"start_date": formatted_start_date,
"end_date": formatted_expiration,
"exp": formatted_expiration,
"ivl": "57600000"
}
response = requests.get(url, headers=headers, params=querystring)
response_json = response.json()
df = wrangle_response(response_json)
df.to_csv(f"C:/Users/anton/Documents/SGH/VI_semestr/RNDs/data/EOD/SPY_option_chain_expiration_{formatted_expiration}.csv")
if not df.empty:
all_data = pd.concat([all_data, df], ignore_index=True)
successful_fetches += 1
print(f"Number of fetches: {successful_fetches}")
all_data.to_csv("C:/Users/anton/VSCode_projects/RNDs/data/SPY_option_chain_expirations3.csv", index=False)
except Exception as e:
print(f"An error occurred: {e}. Data saved up to the last successful fetch.")
break
print(f"Data fetching completed with {successful_fetches} successful fetches.")
return all_databulk_quotes = get_bulk_at_time_quotes("2020-03-20", root = "SPY")And voila. We have it and can proceed further to check it and see how it looks. When doing your analysis, remember to look for the monotonicity of put and call prices. It is a proxy for liquid and illiquid contracts, and the latter might introduce unnecessary noise.
The other factor is standard expirations. There are a lot of expirations that are non-standard, and the options in these expirations could be not liquid enough. It might be a good idea to exclude those dates from the analysis.
Also, check and eliminate any butterfly or calendar spread arbitrage.
Here is the code for further cleaning:
def ms_to_datetime(milliseconds):
seconds, milliseconds = divmod(milliseconds, 1000)
minutes, seconds = divmod(seconds, 60)
hours, minutes = divmod(minutes, 60)
return dt.time(hour=hours, minute=minutes, second=seconds)
def wrangle(df):
df = df.copy()
df["hour"] = df.apply(lambda row: ms_to_datetime(row["ms_of_day"]), axis=1)
df["expiration"] = pd.to_datetime(df["expiration"], format='%Y%m%d')
df["date"] = pd.to_datetime(df["date"], format='%Y%m%d')
df["datetime"] = df.apply(lambda row: pd.Timestamp.combine(row["date"], row["hour"]), axis=1)
df["strike"] = df["strike"]/1000
df["ms_of_day"] = df["ms_of_day"]/ (60*60*1000)
df["mid"] = (df["bid"] + df["ask"]) / 2
df["DTE"] = (df['expiration'] - df['date']).dt.days#
df["tau"] = df["DTE"] / 365
df["maturity_days"] =(df['expiration'] - df['date']).dt.days
columns_to_drop = ["bid_size", "bid_exchange", "bid_condition", "ask_size", "ask_exchange", "ask_condition", "ms_of_day"]
df.drop(columns = columns_to_drop, inplace = True)
return df
#Eliminating all quotes with 0.01 USD bid/ask
def adjust_for_modeling(df):
df = df[df["ask"] > 0.01]
df = df[df["bid"] > 0.01]
df["mid"] = round(df["mid"], 2)
df = df[df["mid"] >= 1/8]
df = df.iloc[:,1:]
return dfFetching SPY prices:
spy_min = min(df_model["date"]).strftime("%Y-%m-%d")
spy_max = max(df_model["date"]).strftime("%Y-%m-%d")
spy = yf.download("SPY", start = spy_min, end = spy_max)
df_spy_adj = spy[["Adj Close", "Close"]]
df_spy_adj.index.name = "date"
df_spy_adj.rename(columns={'Adj Close': 'adj_close', "Close": "close"}, inplace=True)
df_close = pd.merge(df_model, df_spy_adj, on='date', how='left')
df_close = df_close.dropna()Interest rate fetching:
from openbb_terminal.sdk import openbb
openbb.fixedincome.fed("daily", start_date= "2019-07-03")
fed_rates.index.name = "date"
fed_rates.rename(columns={'DFF': 'rate'}, inplace=True)
df_rate = pd.merge(df_close, fed_rates, on='date', how='left')Transformation of the options prices with implied futures prices
Now we are doing an operation, which will help us preserve as much information as possible. Sometimes, people choose between puts and calls. However, there are ITM and OTM regions. We want OTM as they are more luqid and have more information. With the implied futures price approach, we get more information than working just with puts or calls. We effectively exclude ITM options by incorporating just OTM ones. These are more frequently traded for many reasons (e.g., hedging, convexity).
With this approach, you have to calculate implied futures prices to get the spot right, then calculate the prices from put-call parity for ITM options. This is the approach presented in the paper by Ait-Sahalia and Andrew Lo (1998).
The first step is to retrieve the future price for each maturity.
For the derivation of implied futures prices from put-call parity, we have to use the most reliable contract, and these are ATM puts and calls. With the derived future price, we can then plug the future price into the put-call parity equation, which will return the prices of ITM options. Consequently, this method also allows us to circumvent three issues: avoid unreliable data (ITM contracts), have the future/spot price at the time of the recorded prices (we avoid the problem of last traded prices), and circumvent the task of forecasting the correct future dividend yield. We have results, which allow us to drop either put or call prices without any loss of information. Here is the code:
def ITM_valuation_by_parity_put(df):
unique_expirations = df["expiration"].unique()
unique_expirations = unique_expirations
all_results = []
cuonter = 0
for exp_uni in unique_expirations:
# Search for ATM strikes
df_exp = df[df["expiration"] == exp_uni]
df_days = df_exp["date"].unique()
for day in df_days:
df_day = df_exp[df_exp["date"] == day]
Spot = df_day["close"].iloc[0]
ir_uniq = df_day["rate"].iloc[0]
dte_uniq = df_day["tau"].iloc[0]
strike_location = (df_day["strike"] - Spot).abs().sort_values().index[0]
ATM_strike = df_day.loc[strike_location]["strike"]
ATM_strike = np.ceil(ATM_strike) if (ATM_strike - np.floor(ATM_strike)> 0.1 ) else ATM_strike
df_puts = df_day[df_day["right"] == "P"]
df_calls = df_day[df_day["right"] == "C"]
liquid_calls = df_calls[(df_calls["strike"] > ATM_strike)]
illiquid_calls = df_calls[(df_calls["strike"] <= ATM_strike)]
liquid_puts = df_puts[(df_puts["strike"] <= ATM_strike)]
illiquid_puts = df_puts[(df_puts["strike"] > ATM_strike)]
if not df_day[(df_day["strike"] == ATM_strike) & (df_day["right"] == "P")].empty and not df_day[(df_day["strike"] == ATM_strike) & (df_day["right"] == "C")].empty:
put_implied = df_day[(df_day["strike"] == ATM_strike) & (df_day["right"] == "P")]["mid"].iloc[0]
call_implied = df_day[(df_day["strike"] == ATM_strike) & (df_day["right"] == "C")]["mid"].iloc[0]
#print(put_implied,call_implied)
implied_future = call_implied - put_implied + ATM_strike * np.exp(-ir_uniq * dte_uniq)
liquid_calls["put_ITM"] = liquid_calls["mid"] + liquid_calls["strike"] * np.exp(-ir_uniq*dte_uniq) - implied_future
liquid_calls["mid"] = liquid_calls["put_ITM"]
liquid_calls.drop(columns = "put_ITM", inplace = True)
liquid_calls["right"] = "P"
liquid_calls[["bid", "ask"]] = None
# for comparison
illiquid_puts_for_merge = illiquid_puts[['strike', 'mid']].copy()
illiquid_puts_for_merge.rename(columns={'mid': 'illiquid_put_mid'}, inplace=True)
liquid_calls = pd.merge(liquid_calls, illiquid_puts_for_merge, on='strike', how='left')
liquid_calls['mid'].fillna('NA', inplace=True)
result_df = pd.concat([liquid_puts, liquid_calls], ignore_index=True)
all_results.append(result_df)
else:
print(f"No ATM strike for put/call at {exp_uni} expiration. Skipping...")
continue
final_result = pd.concat(all_results, ignore_index=True)
return final_result
Construction of Implied Volatility Surface and interpolation
Now, when we have all the prices, we can proceed with volatility surface interpolation. I won’t delve into the details about what implied volatility is, but in simple words, implied volatility is the volatility input in the Black-Scholes equation, which yields the market price. If you are not familiar with this concept, check on the internet; there are tons of resources.
There are dozens of volatility models, such as local volatility, stochastic ones, and others. I don’t have the quotes to the right of the spot, so we have to choose between not paramatrizing RND and volatilites (we wouldn’t get the right RND tail then) or paramatrizing one of them. I am more comfortable with interpolating in the volatility space, so I will use the SABR model. The SABR stands for “stochastic alpha, beta, and rho,” which are parameters used in the model. It is a model used mostly for interest rate derivatives, but it is correct to use it in this case as well. Note that volatility models differ depending on the asset classes, e.g., in FX markets Vanna-Volga approach is used the most, and in equity markets, local volatility and CEV models are popular.
Nevertheless, I will use the SABR model due to its reliability and wide-spread use by practitioners. It yielded wonderful fits and RND values (more about that later).
from py_vollib.black_scholes_merton.implied_volatility import implied_volatility
def imp_vol_calc_P(row):
S_0 = row["close"]
K = row['strike']
r = row["rate"]
q = 0.0
tau = row['tau']
mid = row['mid']
try:
iv = implied_volatility(mid,S_0, K, tau, r, q, "p")
return iv
except Exception as e:
# Return a NaN or some error indication if calculation fails
return np.nan
df_parity["IV"] = df_parity.apply(imp_vol_calc_P, axis=1)
df_parity["IV"].isna().sum()/len(df_parity["IV"])
# 0.001714122311542935 (Very good!)
from pysabr import Hagan2002LognormalSABR
from pysabr import hagan_2002_lognormal_sabr
def calibrate_sabr(dataframe, beta=0.5):
results = {}
expiries = dataframe["expiration"].unique()
for expiry in expiries:
results[expiry] = {}
sub_df = dataframe[dataframe["expiration"] == expiry]
dates_unique = sub_df["date"].unique()
for date in dates_unique:
sub_df_date = sub_df[sub_df["date"] == date]
sub_df_date = sub_df_date.sort_values(by=["strike"])
strikes = sub_df_date["strike"].values
vols = sub_df_date["IV"].values
close = sub_df_date["close"].iloc[0]
rate = sub_df_date["rate"].iloc[0]
tau = sub_df_date["tau"].iloc[0]
root = sub_df_date["root"].iloc[0]
right = sub_df_date["right"].iloc[0]
datetime = sub_df_date["datetime"].iloc[0]
maturity_days = sub_df_date["maturity_days"].iloc[0]
moneyness = sub_df_date["moneyness"].iloc[0]
forward = close * np.exp(rate*tau)
if len(strikes) > 0 and np.all(vols):
sabr_model = Hagan2002LognormalSABR(f=forward, t = tau, beta=beta)
alpha, rho, volvol = sabr_model.fit(strikes, vols)
results[expiry][date] = {
'alpha': alpha,
'rho': rho,
'volvol': volvol,
'beta': beta,
'expiration': expiry,
'date': date,
'close': close,
'rate': rate,
'tau': tau,
'root': root,
'right': right,
'datetime': datetime,
'maturity_days': maturity_days,
'moneyness': moneyness
}
return results
sabr_calibrated = calibrate_sabr(df_sabr_input)
def extrapolate_sabr_smile(params_dict, strikes_to_left=150, strikes_to_right=150):
all_results = []
for expiry, dates_params in params_dict.items():
for date, params in dates_params.items():
close = params['close']
beta = params.get('beta', 0.5)
rate = params['rate']
tau = params['tau']
strikes = np.arange(round(close - strikes_to_left, 1), round(close + strikes_to_right, 1), 0.1)
forward = close * np.exp(rate * tau)
calibrated_vols = [
hagan_2002_lognormal_sabr.lognormal_vol(k, forward, tau, params['alpha'], beta, params['rho'], params['volvol']) * 100
for k in strikes
]
extended_params = [{**params, 'strike': k, 'IV': vol} for k, vol in zip(strikes, calibrated_vols)]
results_df = pd.DataFrame(extended_params)
results_df.drop(columns=["alpha", "rho", "beta", "volvol"], inplace=True)
all_results.append(results_df)
final_results_df = pd.concat(all_results, ignore_index=True)
return final_results_df
df_sabr_extrapolated = extrapolate_sabr_smile(sabr_calibrated)The perfection of this volatility smile makes me smile as well.
Interpolating with splinces for comparison:
Conversion of implied volatilities back to prices and implementation of the Breeden-Litzenberger formula
Having interpolated the entire volatility surface and having granulated curves for each maturity, this step involves applying the Breeden-Litzenberger formula to gain insights into the market’s expectations of future price movements by calculating RNDs.
By reverse-engineering the Black-Scholes model, one can gain valuable insights into implied distributions and market participants’ expectations, which can be highly informative for investors, traders, or risk professionals. It allows more detailed analysis and insight into market participants’ expectations, market sentiment, or option pricing through the integration of these densities.
To obtain these densities, we have to differentiate the prices twice. The essential formulas are below:
After performing the calculations on market prices, the values obtained are our RND values. As we can see, they look just like normal probability density functions. With expiration, the tails are getting fatter, and the center of the RNDs is shifting toward higher values.
The area below the densities is 1 +/- 0.0001, which is extremely good. I was literally surprised by the accuracy. And they are pretty too.
And the code:
def wrangle_for_interpolation(df):
df = df[["date", "root","expiration",
"strike", "right", "datetime",
"maturity_days", "tau",
"rate", "IV", "close"]]
df["close"] = round(df["close"], 2)
return df
def BS_d1(S, K, r, q, sigma, tau):
return (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * tau) / (sigma * np.sqrt(tau))
def BS_d2(d1, sigma, tau):
return d1 - sigma * np.sqrt(tau)
def BS_price(S, K, r, q, sigma, tau):
d1 = BS_d1(S, K, r, q, sigma, tau)
d2 = BS_d2(d1, sigma, tau)
return K * np.exp(-r * tau) * norm.cdf(-d2) - S * np.exp(-q * tau) * norm.cdf(-d1)
df_sabr_wrangled['bs_price'] = BS_price(df_sabr_wrangled['close'],
df_sabr_wrangled['strike'],
df_sabr_wrangled['rate'], 0,
df_sabr_wrangled['IV'],
df_sabr_wrangled['tau'])
def create_rnd(df):
df = df.sort_values(by = "strike")
prices = df["bs_price"].values
strikes = df["strike"].values
tau = df["tau"].values
r = df["rate"].values
first_grad = np.gradient(prices, strikes)
second_grad = np.gradient(first_grad, strikes)
return strikes, second_grad * np.exp(tau * r)
def create_rnd_multiple(df):
unique_exp = df["expiration"].unique()
result_dfs = []
for exp in unique_exp:
df_exp = df[df["expiration"] == exp]
df_exp_dates = df_exp["date"].unique()
for date in df_exp_dates:
df_day = df_exp[df_exp["date"] == date]
close =df_day["close"].iloc[0]
tau = df["tau"].iloc[0]
r = df["rate"].iloc[0]
rnd_strikes, rnd_values = create_rnd(df_day)
result_dfs.append({
"strike": rnd_strikes,
"rnd": rnd_values,
'expiration': exp,
'date': date,
"close": close,
'tau': tau,
"rate": r,
})
df_final = pd.DataFrame(result_dfs)
df_final = df_final.explode(['strike', 'rnd']).reset_index(drop=True)
return df_final
def area_under_rnds(df):
unique_exp = df["expiration"].unique()
results = []
for expiration in unique_exp:
df_exp = df[df["expiration"] == expiration]
df_day_unique = df_exp["date"].unique()
for day in df_day_unique:
df_day_exp = df_exp[df_exp["date"] == day]
interval_width = df_day_exp['strike'].diff().fillna(0)
area_under_curve = round((df_day_exp['rnd'] * interval_width).sum(), 6)
results.append({'Expiration': expiration, "date": day, 'Area Under Curve': area_under_curve})
result_df = pd.DataFrame(results)
return result_dfInterpretation
The risk-neutral densities are indeed risk-neutral, which means that these are not real-world probabilities. That means that these probabilities are correct in a world indifferent to risk. This is a correct framework for pricing derivatives; however, a straightforward translation to our world could be incorrect. If you wanted to translate it to the real world, you would have to estimate the investors’ utility function. It is another step that requires assumptions, which, in the end, could be unreliable. Too many assumptions are problematic, they can introduce unnecessary problems.
Nevertheless, there is a link. Investors’ sentiments are reflected in prices, which, converted to densities, should be interpretable. Consequently, changes in RNDs are an arbitrage-free reflection of their preferences and the relative weight attached to the upside and downside. Their “skin in the game,” measured by the prices, can somehow be attached to these probabilities.
There are also some applications of RNDs in terms of risk management purposes or trading. We can trade kurtosis and skewness with them (sounds fancy, but this boils down to agreeing or not with their shape and placing appropriate trades). We can think that, for example, the probability density function looks too narrow and thin-tailed. We can bet on it. To extract the 2nd and 4th moments (mean and kurtosis) of such distribution, we can sell ATM straddles and buy OTM straddles in greater quantity, or, for example, buy calendar spreads. Both trades allow us to bet on the price distribution.
One can also compute the kurtosis or skewness time series and use it to generate signals. Too small a kurtosis at some point might mean that investors are anticipating little changes in prices and no big move. If you are contrarian, you can bet on it.
Note that RNDs are also histograms, which means that they do not show path- dependence. This implies that if you are a dynamic hedger, path dependence is important to you. Hedging the delta introduces that, and RND can be useless as it takes into consideration terminal prices and not the relation between volatility and price level.
Conclusion
I am content with the present RND values, given the primary objective — extracting the densities. The RNDs are consitent and reliable. We can do loads of things with them and I invite to experimetning with the code yourself.
I hope you enjoyed this article. If you have any questions, you can reach out to me on LinkedIn.
