Why is my Risk Neutral Density recovery failing?

Question

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}')

Answer 1

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