In [3]:
plot_return_distribution_with_bernoulli_events(
background_std = .1,
event_prob = 0.0,
event_std = .2,
background_mean = 0.0,
event_mean = 0.0
)
I spent my junior summer at a FICC Options desk in Chicago. There, one of the most unintuitive aspects of pricing options was skew: the observation that implied volatility of an asset will vary across different strike prices with the same time to expiry.
The most basic pricing model for options is Black Scholes, which assumes that asset returns are normally distributed with a known mean and standard deviation. Real life is not that neat -- it is intuitive that an asset's returns will be more volatile on certain days (e.g. earnings, announcements, during corporate actions) than others. The implied volatility of an option reflects the market's expectation for how much volatility the asset will incur over the life of an option. Thus, for any given strike price, different expiries for the same asset will have different implied volatilities.
Still, I struggled to understand how over the exact same duration of time, different strikes could have different implied volatilities in a bizarre superposition of standard deviations. I heard explanations involving how far OTM puts were in demand from hedgers, and far OTM calls were in demand from people seeking cheap leveraged exposure to assets. This didn't feel like the full picture -- the structure of skew varies significantly between assets, and many assets have fairly stable skew structures across time. If the phenomenon arose from supply and demand, I'd expect the extent of skew to fluctuate a lot more. The explanation that finally satisfied me was that far OTM options would only have value if there was a large move in asset prices. Conditional on there being a large move, the confidence interval over which asset prices could fall also becomes larger.
Here, I want to explore a simple model. I'll only examine options over a single TTE (for simplicity, we can call this a week). Every time period, asset returns will move according to a baseline normal distribution. Additionally, during some of these time periods, there can be an event. The returns of an event will have their own mean and standard deviation, and will be added to the baseline returns. I want to bolster my own intuition about how skew can arise from this simple, plausible setup.
The first section will cover the shape of the distribution of returns under such a model. Ths second section will visualize put skew as these parameters change. The third section will calculate MLE fits of different stocks under this model. In addition to static graphs, I'll add some code to interact with the model by toggling the parameters in a Jupyter notebook for independent exploration.
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm, bernoulli, t, laplace, kstest
from scipy.optimize import brentq
from ipywidgets import interact, FloatSlider, FloatLogSlider, interactive
from typing import Optional, Union, Callable, Tuple
import warnings
import matplotlib.pyplot as plt
def plot_return_distribution_with_bernoulli_events(
background_std: float,
event_prob: float,
event_std: float,
background_mean: float = 0.0,
event_mean: float = 0.0
):
"""
Interactive widget to plot a histogram of returns if daily returns are governed by the following process:
- Returns have a standard background mean/standard deviation
- There is some probability of an event, governed by a Bernoulli random variable
- If there is an event, that is governed by a separate event mean and standard deviation
- Daily returns are N(background_mean, background_std) + event_indicator * N(event_mean, event_std)
This widget plots a histogram of daily returns drawn from Monte Carlo simulation with overlayed MLE fits and their respective Kolmogorov-Smirnov (KS) statistics:
- the MLE normal distribution describing returns
- the MLE Laplace distribution describing returns
Params:
background_std (float) - the unavoidable daily standard deviation in returns
event_prob (float) - the probability of an event happening, on the domain [0, 1]
event_std (float) - given that an event happens, the standard deviation of its impact on returns
background_mean (float) - the mean of background returns
event_mean (float) - the mean of event returns
"""
# Simulate data
background_returns = np.random.normal(background_mean, background_std, 100000)
event_indicator = bernoulli.rvs(event_prob, size=100000)
event_returns = np.random.normal(event_mean, event_std, 100000)
result = background_returns + event_indicator * event_returns
# MLE Normal fit
mu, std = norm.fit(result)
# MLE Laplace fit
laplace_loc = np.mean(result)
laplace_scale = np.std(result) / np.sqrt(2)
# Calculate the Bayes (true) CDF using the known process
background_cdf = norm.cdf(result, background_mean, background_std)
event_cdf = norm.cdf(result, event_mean, event_std)
true_cdf = background_cdf + event_prob * event_cdf
# Calculate the Normal and Laplace CDFs
normal_cdf = norm.cdf(result, mu, std)
laplace_cdf = laplace.cdf(result, loc=laplace_loc, scale=laplace_scale)
# Calculate Kolmogorov-Smirnov statistics
ks_normal = np.max(np.abs(true_cdf - normal_cdf))
ks_laplace = np.max(np.abs(true_cdf - laplace_cdf))
# Plotting
plt.figure(figsize=(10, 6))
# Plot histogram of the results
sns.histplot(result, bins=200, kde=False, stat='density', color='lightblue', label="Returns")
# Plot Normal and Laplace fits
x = np.linspace(min(result), max(result), 1000)
plt.plot(x, norm.pdf(x, mu, std), label=f'Normal (KS: {ks_normal:.2f})', color='blue')
plt.plot(x, laplace.pdf(x, loc=laplace_loc, scale=laplace_scale),
label=f'Laplace (KS: {ks_laplace:.2f})', color='red')
# Add title and labels
plt.legend(title="KS Statistics")
plt.title("Returns Sampled Distribution with KS Statistics")
plt.xlabel("Returns")
plt.ylabel("Density")
plt.grid(False)
plt.show()
plot_return_distribution_with_bernoulli_events(
background_std = .1,
event_prob = 0.0,
event_std = .2,
background_mean = 0.0,
event_mean = 0.0
)
plot_return_distribution_with_bernoulli_events(
background_std = .1,
event_prob = 0.2,
event_std = .2,
background_mean = 0.0,
event_mean = 0.0
)
plot_return_distribution_with_bernoulli_events(
background_std = .1,
event_prob = 0.2,
event_std = .4,
background_mean = 0.0,
event_mean = 0.0
)
plot_return_distribution_with_bernoulli_events(
background_std = .1,
event_prob = 0.44,
event_std = .2,
background_mean = 0.0,
event_mean = 0.0
)
interactive_plot = interactive(
plot_return_distribution_with_bernoulli_events,
background_mean=FloatSlider(min=-0.1, max=0.1, step=0.01, value=0.0, description="Bg Mean"),
background_std=FloatSlider(min=0.0, max=0.5, step=0.01, value=0.1, description="Bg Std"),
event_mean=FloatSlider(min=-0.1, max=0.1, step=0.01, value=0.0, description="Event Mean"),
event_std=FloatSlider(min=0.0, max=0.5, step=0.01, value=0.1, description="Event Std"),
event_prob=FloatSlider(min=0.0, max=1.0, step=0.01, value=0.1, description="Event Prob")
)
interactive_plot
interactive(children=(FloatSlider(value=0.1, description='Bg Std', max=0.5, step=0.01), FloatSlider(value=0.1,…
def price_option(
S: float,
K: float,
baseline_mean: float,
baseline_std: float,
event_mean: float,
event_std: float,
p_event: float,
option_type: str = 'call'
) -> float:
"""
Price an option using a mixed normal distribution model.
Args:
S: Current stock price
K: Strike price
baseline_mean: Mean of baseline normal distribution
baseline_std: Standard deviation of baseline distribution
event_mean: Additional mean if event occurs
event_std: Additional std if event occurs
p_event: Probability of event occurring
option_type: 'call' or 'put'
Returns:
Option price
"""
def _normal_price(mean: float, std: float) -> float:
# Guard against division by zero
if std == 0:
# For zero volatility, option value is max(S*exp(mean) - K, 0) for calls
# or max(K - S*exp(mean), 0) for puts
forward = S * np.exp(mean)
if option_type.lower() == 'call':
return max(forward - K, 0)
else:
return max(K - forward, 0)
d1 = (np.log(S/K) + mean + 0.5 * std**2) / std
d2 = d1 - std
if option_type.lower() == 'call':
return S * np.exp(mean) * norm.cdf(d1) - K * norm.cdf(d2)
else:
return K * norm.cdf(-d2) - S * np.exp(mean) * norm.cdf(-d1)
# When p_event is 0, we only need the baseline case
if p_event == 0:
return _normal_price(baseline_mean, baseline_std)
# When p_event is 1, we only need the event case
if p_event == 1:
combined_mean = baseline_mean + event_mean
combined_std = np.sqrt(baseline_std**2 + event_std**2)
return _normal_price(combined_mean, combined_std)
# Regular case: weighted average of both scenarios
no_event_price = _normal_price(baseline_mean, baseline_std)
combined_mean = baseline_mean + event_mean
combined_std = np.sqrt(baseline_std**2 + event_std**2)
event_price = _normal_price(combined_mean, combined_std)
return (1 - p_event) * no_event_price + p_event * event_price
def calc_implied_vol(
price: float,
S: float,
K: float,
T: float = 1.0,
option_type: str = 'call',
r: float = 0.0,
initial_guess: float = 0.3,
tolerance: float = 1e-5
) -> Optional[float]:
"""
Calculate implied volatility using Brent's method with robust error handling.
Args:
price: Market price of option
S: Current stock price
K: Strike price
T: Time to expiration in years
option_type: 'call' or 'put'
r: Risk-free rate
initial_guess: Initial volatility guess
tolerance: Tolerance for convergence
Returns:
Implied volatility or None if no solution found
"""
forward = S * np.exp(r * T)
# Handle edge cases for zero or negative option prices
if price <= 0:
return 0.0
# Calculate intrinsic value
if option_type.lower() == 'call':
intrinsic = max(0, forward - K)
else:
intrinsic = max(0, K - forward)
# If price is at intrinsic value, return 0 volatility
if abs(price - intrinsic) < tolerance:
return 0.0
def bs_price(sigma: float) -> float:
if sigma < tolerance:
return intrinsic
d1 = (np.log(S/K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
d2 = d1 - sigma * np.sqrt(T)
if option_type.lower() == 'call':
return S * np.exp(-r * T) * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
else:
return K * np.exp(-r * T) * norm.cdf(-d2) - S * np.exp(-r * T) * norm.cdf(-d1)
try:
# Try increasingly wider ranges
ranges = [(0.0001, 1.0), (0.0001, 2.0), (0.00001, 3.0)]
for low, high in ranges:
try:
return brentq(lambda x: bs_price(x) - price, low, high, rtol=tolerance)
except ValueError:
continue
return None
except Exception as e:
warnings.warn(f"Failed to compute implied vol: {str(e)}")
return None
def generate_vol_surface(
pricing_func: Callable,
S: float,
params: dict,
min_strike_mult: float = 0.5,
max_strike_mult: float = 1.5,
num_strikes: int = 300,
clean_nans: bool = True
) -> Tuple[np.ndarray, np.ndarray]:
"""
Generate implied volatility surface with robust error handling.
Args:
pricing_func: Function that prices options
S: Current stock price
params: Dictionary of parameters for pricing_func
min_strike_mult: Minimum strike price as multiple of S
max_strike_mult: Maximum strike price as multiple of S
num_strikes: Number of strike prices to evaluate
clean_nans: If True, interpolate over NaN values
Returns:
Tuple of (strikes, implied_vols)
"""
strikes = np.linspace(S * min_strike_mult, S * max_strike_mult, num_strikes)
# calculate the forward price
baseline_mean, baseline_std, event_mean, event_std, p_event = params.values()
forward = S * np.exp(baseline_mean) * (1-p_event) + S * np.exp(baseline_mean + event_mean) * p_event
implied_vols = []
option_type = "put"
for K in strikes:
try:
price = pricing_func(S, K, option_type=option_type, **params)
# Handle zero or negative prices
if price <= 0:
implied_vols.append(0.0)
continue
vol = calc_implied_vol(price, S, K, option_type=option_type)
implied_vols.append(vol if vol is not None else np.nan)
except Exception as e:
warnings.warn(f"Error at strike {K}: {str(e)}")
implied_vols.append(np.nan)
vols = np.array(implied_vols)
# Clean up NaN values through interpolation if requested
if clean_nans and np.any(np.isnan(vols)):
valid_mask = ~np.isnan(vols)
if np.any(valid_mask):
x = np.arange(len(vols))
valid_indices = x[valid_mask]
valid_vols = vols[valid_mask]
vols = np.interp(x, valid_indices, valid_vols)
return strikes, vols
def plot_vol_smile(
S: float,
baseline_mean: float,
baseline_std: float,
event_mean: float,
event_std: float,
p_event: float
) -> None:
"""
Create plot of implied volatility smile with robust axis scaling.
Args:
S: Current stock price
baseline_mean: Mean of baseline distribution
baseline_std: Std dev of baseline distribution
event_mean: Additional mean if event occurs
event_std: Additional std if event occurs
p_event: Probability of event occurring
"""
params = {
'baseline_mean': baseline_mean,
'baseline_std': baseline_std,
'event_mean': event_mean,
'event_std': event_std,
'p_event': p_event
}
try:
strikes, vols = generate_vol_surface(price_option, S, params)
# Find ATM strike for reference line
forward_strike = S * np.exp(baseline_mean) * (1-p_event) + S * np.exp(baseline_mean + event_mean) * p_event
# Clean any remaining NaN values for max calculation
clean_vols = vols[~np.isnan(vols)]
if len(clean_vols) == 0:
raise ValueError("No valid volatility values found")
# Calculate axis limits
max_vol = np.max(clean_vols)
min_moneyness = min(strikes) / S
max_moneyness = max(strikes) / S
plt.figure(figsize=(10, 6))
# Plot implied volatility line
plt.plot(strikes/S, vols, 'b-', label='Implied Volatility', linewidth=2)
# Add FW reference line
plt.axvline(x=forward_strike/S, color='g', linestyle=':',
label=f'Forward (K={forward_strike:.1f})', alpha=0.5)
# Add baseline std as reference line if it's a pure baseline case
if p_event == 0:
plt.axhline(y=baseline_std, color='r', linestyle='--',
label=f'Baseline σ={baseline_std:.2f}', linewidth=1.5)
# Set axis limits and labels
plt.xlim(min_moneyness, max_moneyness)
plt.ylim(0, max(max_vol * 1.5, baseline_std * 1.5))
plt.xlabel('Moneyness (K/S)')
plt.ylabel('Implied Volatility')
plt.title(f'Implied Volatility Smile (p_event={p_event:.2f})')
# Enhance grid and legend
plt.grid(True, alpha=0.3)
plt.legend(loc='best')
# Add horizontal grid lines at round numbers
max_y = plt.ylim()[1]
grid_steps = np.arange(0, max_y + 0.1, 0.1)
plt.hlines(grid_steps, min_moneyness, max_moneyness,
colors='gray', linestyles=':', alpha=0.2)
plt.show()
except Exception as e:
plt.close() # Clean up any partial plot
warnings.warn(f"Failed to generate plot: {str(e)}")
raise
def create_interactive_plot() -> None:
"""
Create interactive widget controls with reasonable bounds and steps.
"""
return interact(
plot_vol_smile,
S=FloatSlider(value=100.0, min=50.0, max=150.0, step=1.0,
description='Stock Price'),
baseline_mean=FloatSlider(value=0.0, min=-0.2, max=0.2, step=0.01,
description='Base Mean'),
baseline_std=FloatSlider(value=0.2, min=0.0, max=0.5, step=0.01,
description='Base Std'),
event_mean=FloatSlider(value=-0.05, min=-0.2, max=0.2, step=0.01,
description='Event Mean'),
event_std=FloatSlider(value=0.3, min=0.0, max=1.0, step=0.01,
description='Event Std'),
p_event=FloatSlider(value=0.1, min=0.0, max=1.0, step=0.01,
description='Event Prob')
)
plot_vol_smile(
S = 100,
baseline_mean = 0.0,
baseline_std = .2,
event_mean = 0.0,
event_std = 0.0,
p_event = 0.1
)
plot_vol_smile(
S = 100,
baseline_mean = 0.0,
baseline_std = .2,
event_mean = 0.0,
event_std = 0.4,
p_event = 0.1
)
plot_vol_smile(
S = 100,
baseline_mean = 0.0,
baseline_std = .2,
event_mean = -0.2,
event_std = 0.4,
p_event = 0.5
)
plot_vol_smile(
S = 100,
baseline_mean = 0.0,
baseline_std = .2,
event_mean = 0,
event_std = 0.4,
p_event = 0.3
)
create_interactive_plot()
interactive(children=(FloatSlider(value=100.0, description='Stock Price', max=150.0, min=50.0, step=1.0), Floa…
<function __main__.plot_vol_smile(S: float, baseline_mean: float, baseline_std: float, event_mean: float, event_std: float, p_event: float) -> None>
As shown above, the existence of a volatile event in excess of a day's baseline volatility implies the existence of skew in the "fair" IV of options at different strike prices.
I don't think this model is useful for justifying whether call skew or put skew should be larger -- there are natural distortions in this simulation because a log-normal distribution cannot go below 0, and the avg price of the asset in time period 1 is affected by the mean return AND the standard deviation of the returns (this is also why phenomena like [Shannon's Demon] (https://portfoliocharts.com/2022/04/12/unexpected-returns-shannons-demon-the-rebalancing-bonus/) arise in portfolio management).
Notice that the skew of the volatility surface peaks when the probability of an event is 50%. After that point, a more probable event just becomes a new base state, and the extent of skew in the tails starts decreasing.
The mathematical intuition behind this skew is that, if the options in the tails become valuable, it is likely that the asset experienced an event that has dramatically widened the cone of possible prices upon expiry. To compensate for this risk, the options in the tail should have higher implied volatility.