0
$\begingroup$

I'm working on a project to recover a known Risk Neutral Density from option prices, using the Breeden-Litzenberger formula (assuming a continuum of option_price(strike_price), the second derivative of an option's price with respect to the strike price is equal to the risk-neutral probability density function of the underlying asset's price).

I'm testing my formula with a uniform RND, but the recovered "density" ends up being 100 times larger than expected (100*densityR). The issue occurs no matter what density I use (uniform or an approximated-normal).

Here's what I'm doing:

  1. I generate a linearly spaced returns vector and calculate StrikePrices.

  2. I generate a uniform discrete density densityR.

  3. I calculate option prices for puts (option price = discounted expected payoff under (strike prices, densityR) )

  4. I calculate the second derivative of option prices with respect to strike prices to recover the density.

Both MATLAB and Python implementations yield similar discrepancies. Here are some screenshots: https://imgur.com/a/s66kOY8

Why is the recovered "density" 100 times bigger than expected? Where am I making a mistake? I would greatly appreciate any assistance.

MATLAB version:

r=0.02;
tau=0.25;
numElements = 10000;
currentPrice=49.75;

% Create returns as a linearly spaced vector from 0 to 2
returns = linspace(0, 2, numElements)';

StrikePrices = currentPrice * returns; 

% Create densityR as a uniform distribution
densityR = ones(numElements, 1) * (1 / numElements);

% Calculate the payoff for each option price
optionPrices = zeros(numElements, 1); % Initialize optionPrices array

% Adjusted loop to calculate option prices
for i = 1:numElements
    allPayoffsForCurrentStrike = max(0, StrikePrices(i) - StrikePrices);
    optionPrices(i) = sum(allPayoffsForCurrentStrike .* densityR) * exp(-r * tau);
end

optionPrices=optionPrices';

% Calculate the first derivative of optionPrices with respect to StrikePrices
dy_dx = gradient(optionPrices, StrikePrices);

% Calculate the second derivative of optionPrices with respect to StrikePrices
d2y_dx2 = gradient(dy_dx, StrikePrices);

first_derivative=dy_dx;
second_derivative=d2y_dx2;
sum_of_derivative=sum(second_derivative)

figure;
hold on; 
plot(dy_dx, 'o-', 'DisplayName', 'First Derivative'); 
plot(second_derivative, 'o-', 'DisplayName', 'Second Derivative'); 
plot(optionPrices, 'o-', 'DisplayName', 'Option Prices'); 
plot(StrikePrices, 'o-', 'DisplayName', 'Strike Prices'); 
plot(densityR, 'o-', 'DisplayName', 'Density R'); 
hold off; 
title('Second Derivative of Option Prices and Density R');
xlabel('Index');
ylabel('Value');
legend show; 
grid on; 

set(gca, 'YScale', 'log'); 

Python version:

import numpy as np
import matplotlib.pyplot as plt

# Given parameters
r = 0.02
tau = 0.25
numElements = 10000
currentPrice = 49.75 

# Create returns as a linearly spaced vector from 0 to 2
returns = np.linspace(0, 2, numElements)

# Calculate StrikePrices
StrikePrices = currentPrice * returns

# Create densityR as a uniform distribution
densityR = np.ones(numElements) * (1 / numElements)

# Calculate the option prices
optionPrices = np.zeros(numElements)
for i in range(numElements):
    allPayoffsForCurrentStrike = np.maximum(0, StrikePrices[i] - StrikePrices)
    optionPrices[i] = np.sum(allPayoffsForCurrentStrike * densityR) * np.exp(-r * tau)

# Calculate the first derivative of optionPrices with respect to StrikePrices
dy_dx = np.gradient(optionPrices, StrikePrices)

# Calculate the second derivative of optionPrices with respect to StrikePrices
d2y_dx2 = np.gradient(dy_dx, StrikePrices);

# Sum of the second derivative
sum_of_derivative = np.sum(d2y_dx2)

# Plotting
plt.figure()
plt.plot(dy_dx, 'o-', label='First Derivative')
plt.plot(d2y_dx2, 'o-', label='Second Derivative')
plt.plot(optionPrices, 'o-', label='Option Prices')
plt.plot(StrikePrices, 'o-', label='Strike Prices')
plt.plot(densityR, 'o-', label='Density R')
plt.title('Second Derivative of Option Prices and Density R')
plt.xlabel('Index')
plt.ylabel('Value')
plt.yscale('log')
plt.legend()
plt.grid(True)
plt.show()

print(f'Sum of the second derivative: {sum_of_derivative}')
$\endgroup$

1 Answer 1

1
$\begingroup$

You need to rescale the density of $S_T$, to recover the density of your return factor $R$. For example, if $R \sim U(0,2)$, then $S_T \sim U(0,2S_0)$.

As you have written it, the variable densityR does not denote the density of $R$ at each point, but rather the probability that $R$ lies in each partition interval. If you wish to compute the option-prices using Riemann sums, you could write

# Create densityR as a uniform distribution
densityR = np.ones(numElements) / 2
densityS = densityR / (currentPrice) # For the uniform distribution
# Compute asset density over StrikePrices

# Calculate the option prices
optionPrices = np.zeros(numElements)
for i in range(numElements):
    allPayoffsForCurrentStrike = np.maximum(0, StrikePrices[i] - StrikePrices)
    optionPrices[i] = np.sum(allPayoffsForCurrentStrike[:-1] * densityS[:-1] * np.diff(StrikePrices)) * np.exp(-r * tau)

If you compare densityS to the estimated density d2ydx2, you will find that they agree except at the very edges of their range (due to numerical errors).

$\endgroup$
12
  • 1
    $\begingroup$ Thank you, simply doing densityS = densityR / (2*(currentPrice)), and then using it instead of densityR in my original option price calculation ensures that (the second derivative of optionPrices with respect to StrikePrices) ≈ densityR. However, densityS sums to ~0.01, not 1. I've never done these sort of discretization problems before, so I might be missing something obvious. Can you tell me, and/or link me to further reading material, regarding why this rescaling is justified and necessary? $\endgroup$
    – v.y.
    Feb 28 at 21:33
  • 1
    $\begingroup$ I don't understand why densityS = densityR / (2*(currentPrice)) gives us the density of ST across the StrikePrices range. How do you derive that from "the probability that R lies in each partition interval = densityR"? $\endgroup$
    – v.y.
    Feb 29 at 20:24
  • 1
    $\begingroup$ I understand that densityS = densityR / (2*(currentPrice)) follows from the change of variables in probability. However, you said "densityR does not denote the density of R at each point, but rather the probability that R lies in each partition interval" - if that is so, how could we apply the change of variables rules, and produce a "discrete sample of the density function" - when densityR is apparently not a density? $\endgroup$
    – v.y.
    Mar 1 at 20:00
  • 1
    $\begingroup$ Also, how would I specify densityS directly? Even if I do StrikePrices = np.linspace(0, currentPrice*2, numElements) and densityS = np.ones(numElements) * (1 / numElements), I still get the same issue as before (imgur.com/a/jPpMvfA). $\endgroup$
    – v.y.
    Mar 1 at 20:03
  • 1
    $\begingroup$ I wrote densityS = scipy.stats.uniform(loc=0, scale=2*currentPrice), and use it as pdf_values=densityS.pdf(StrikePrices) ... probability_of_lying_in_each_interval=np.insert(pdf_values[1:]*np.diff(StrikePrices),0,0) ... optionPrices[i] = np.sum(allPayoffsForCurrentStrike * probability_of_lying_in_each_interval) * np.exp(-r * tau). When I do this, np.sum(probability_of_lying_in_each_interval)=1, and I do manage to recover the original PDF (imgur.com/a/1ty0MLm), and also the probability_of_lying_in_each_interval. Thanks. $\endgroup$
    – v.y.
    Mar 4 at 21:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.