# Why is my Risk Neutral Density recovery failing?

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

% Calculate the second derivative of optionPrices with respect to 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

# Calculate the second derivative of optionPrices with respect to 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}')


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).

• 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?
– v.y.
Commented Feb 28 at 21:33
• 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"?
– v.y.
Commented Feb 29 at 20:24
• 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?
– v.y.
Commented Mar 1 at 20:00
• 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).
– v.y.
Commented Mar 1 at 20:03
• 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.
– v.y.
Commented Mar 4 at 21:37