# Estimation of Risk-Neutral Densities Using Positive Convolution Approximation - Python

I'm trying to estimate the risk-neutral density through positive convolution approximation (introduced by Bondarenko 2002: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=375781). I'm currently struggling to implement a computational algorithm in Python solving the following optimization problem:

$$\hat{f}(u) := \sum_{j}a_j \phi(u - z_j)$$

$$\min_{a}\sum_{i=1}^n\left(P_i - \int_{-\infty}^{x_i}\left(\int_{-\infty}^y \hat{f}(u) du\right)dy\right)^2 \\ s.t. \sum_{j}a_j=1\\ \forall a_j\geq0$$

where: $$P_i$$ is the observed put price with strike $$x_i$$, $$\phi$$ refer to a rescaled standard normal distribution and $$z_j$$ represents points on an equally-spaced grid between the smallest and largest observed strikes $$x$$.

It is straightforward that one should use a standard quadratic program to solve the problem. However, I don't know how to handle the fact that the $$a$$'s are inside a function of $$u$$, which itself is inside a double integral.

Does anyone already implemented positive convolution approximation to estimate the risk-neutral density in Python?

Otherwise, could someone show me how to code an optimization problem with a double integral, such as for example:

$$\min_a\int_{-\infty}^{x}\left(\int_{-\infty}^y \hat{f}(u) du\right)dy \\ \hat{f}(u) := \sum_{j}a_j (u - z_j)^2$$

Thanks for the help!

EDIT

Update: Thanks to the comments and the answer of Attack68. I was able to implement the following code:

import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm

# Compute f_hat
def f(u, y, *args):
a = args[0]
z = args[1]
h = args[2]
j = len(a)
#    print(np.sum(a * norm.pdf(np.tile(u, [j,1]).transpose(), z, h), axis=1))
return np.sum(a * norm.pdf(np.tile(u, [j,1]).transpose(), z, h), axis=1)

# Compute double integral
def DI(a, b, z, h):
#    print(dblquad(f, -10000, b, lambda x: -10000, lambda x: x, args=(a, z, h))[0])
return dblquad(f, -np.inf, b, lambda x: -np.inf, lambda x: x, args=(a, z, h))[0]

def sum_squared_pricing_diff(a, P, X, z, h):
total = 0
for i in range(0, len(P)):
p = P[i]
x = X[i]
total += (p - DI(a, x, z, h)) ** 2

# P is an array of vector put option prices
P = [0.249999283, 0.43750315, 0.572923413, 0.760408034, 0.94790493, 1.14584317,
1.458335038, 1.77083305, 2.624999786, 3.812499791, 5.377596753, 8.06065865,
10.74376984, 14.88873497, 19.88822895]

# X is the vector of the corresponding strikes of the put options
X = [560, 570, 575, 580, 585, 590, 595, 600, 605, 610, 615, 620, 625, 630, 635]

# z is the equally-spaced grid
z = np.linspace(0, 1000, 20)
# h arbitrarily chosen
h = 0.5
# initial guess of a
a_0 = np.ones(len(z)) / len(z)

constraints = ({'type': 'eq', 'fun': lambda a: 1 - np.sum(a)},)
bounds = (((0,None),)*len(z))

sol = minimize(sum_squared_pricing_diff, a_0, args=(P, X, z, h), method='SLSQP', constraints=constraints, bounds=bounds)
print(sol)


which returns the following warning and has difficulty to converge:

IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties.  If the position of a
local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges.  Perhaps a special-purpose integrator should be used.
warnings.warn(msg, IntegrationWarning)


Following a stack overflow post I will try to use nquad instead of dblquad. I will post further progress.

EDIT 2 Update: Using the insights from Attack68's second answer, I was able to estimate the RND in an "efficient" way (probably it can be further improved):

import pandas as pd
import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm
import matplotlib.pyplot as plt
import math

###############################################################################
# Define required functions to describe the optimization problem
###############################################################################

# Double integral transformed
def sum_j_aK(a, x, z, h):
j = len(a)
loc = z
scale = h
x_normalized = (np.ones(j)*x - loc) / scale
K_j = (x_normalized*norm.cdf(x_normalized) + np.exp(-0.5*x_normalized**2)/((2*np.pi)**0.5)) * scale
return np.sum(a*K_j)

# Minimization problem
def sum_squared_pricing_diff(a, P, X, z, h):
total = 0
for i in range(0, len(P)):
p = P[i]
x = X[i]
total += abs(p - sum_j_aK(a, x, z, h))

###############################################################################
# Input required to solve the optimization problem
###############################################################################

# P is an array of vector put option prices
P = [0.249999283, 0.43750315, 0.572923413, 0.760408034, 0.94790493, 1.14584317,
1.458335038, 1.77083305, 2.624999786, 3.812499791, 5.377596753, 8.06065865,
10.74376984, 14.88873497, 19.88822895]

# X is the vector of the corresponding strikes of the put options
X = [560, 570, 575, 580, 585, 590, 595, 600, 605, 610, 615, 620, 625, 630, 635]

# h and j can be chosen arbitrarily
h = 4 # the higher h the smoother the estimated risk-neutral density
j = 50 # the higher the slower the optimization process

###############################################################################
# Solving the optimization problem
###############################################################################

# z is the equally-spaced grid
z = np.linspace((int(math.floor(min(X) / 100.0)) * 100), (int(math.ceil(max(X) / 100.0)) * 100), num=j)

# initial guess of a
a_0 = np.ones(j) / j

# The a vector has to sum up to 1
constraints = ({'type': 'eq', 'fun': lambda a: 1 - np.sum(a)},)

# Each a has to be larger or equal than 0
bounds = (((0,None),)*j)

sol = minimize(sum_squared_pricing_diff, a_0, args=(P, X, z, h), method='SLSQP', constraints=constraints, bounds=bounds)
print(sol)

###############################################################################
# Visualize obtained risk-neutral density (rnd)
###############################################################################
n = 500
a_sol = sol.x
s = np.linspace(min(X)*0.8, max(X)*1.2, num=n)

rnd = pd.DataFrame(np.sum(a_sol * norm.pdf(np.tile(s, [len(a_sol),1]).transpose(), z, h), axis=1))
rnd.index = s

plt.figure()
plt.plot(rnd)

• So the $a_j$ terms are just constants? If yes, then just pull them out of the double integral? – LocalVolatility Feb 19 at 17:02
• @LocalVolatility Thanks. Yes, they are just constants. Indeed, could be a way to tackle the problem. However, this would imply that I would need to vectorize the integrals over an array and therefore could not use scipy.integrate.dblquad since it employs an adaptive algorithm. Thus, I would need to use np.trapz, scipy.integrate.simps, or scipy.integrate.romb and not be allowed to use a function object. Right? – William Burknecht Feb 19 at 21:10

Looking at this a second time, your implementation is probably not very inefficient since for; $$\quad \hat{f}(u) := \sum_{j}a_j \phi(u - z_j), \quad \int_{-\infty}^y \hat{f}(u) du = \sum_j a_j \Phi(y)\;,$$

where $$\Phi(y)$$ is the transformed cumulative normal distribution function of $$y$$. This function already exists as scipy.stats.norm.cdf and is optimised.

Then you have; $$\int_{-\infty}^{x_i} \sum_j a_j \Phi(y) dy = \sum_j a_j \int_{-\infty}^{x_i} \Phi(y) dy = \sum_j a_j K_j \;.$$

These $$K_j$$ are just constants from the point of view of the $$a_j$$ optimisation so can be precomputed, and since you have scipy.stats.norm.cdf it seems you only need to use quad and not dblquad. You might also be interested to know that for a standard normal cdf (see https://math.stackexchange.com/questions/2040980/solving-approximating-integral-of-standard-normal-cdf ):

$$\int_{-\infty}^{x} \Phi(y) dy = x \Phi(x) + \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} \;.$$

$$\min_{a_j} \sum_i \left ( P_i - \sum_j a_j K_j \right )^2 \; s.t. constraints$$

With a bit more analysis I wonder if that can't be further simplified...

scipy.integrate.dblquad solves the following iterated integral problem:

$$I(\alpha,\beta) = \int_{\alpha}^{\beta} \int_{g(x)}^{h(x)} f(x,y)\; dy \;dx$$

If I re-express your variables to be consistent with the scipy documentation, and suppose I changed your integral to be finite:

$$\min_\alpha \int_{-\infty}^{\beta}\left(\int_{-\infty}^x \hat{f}(y) dy\right)dx \\ \hat{f}(y) := \sum_{j}a_j e^{(y - z_j)}$$

Then would it not be as simple as:

import numpy as np
from scipy.optimize import minimize

def f(y, x, *args):
a = args[0]
z = args[1]
return np.sum(a * np.exp(y - z))

def I(a, b, z):
return dblquad(f, -np.inf, b, lambda x: -np.inf, lambda x: x, args=(a, z))[0]


And to engage the optimisation you might try:

n = 5; b = 2; z = np.array([-2,0,3,0.5,2.2])
a_0 = np.ones(n) / n

constraints = ({'type': 'eq', 'fun': lambda a: 1 - np.sum(a)},)
bounds = (((0,None),)*n)

sol = minimize(I, a_0, args=(b, z), method='SLSQP', constraints=constraints, bounds=bounds)
print(sol)

>>>    fun: 0.3678794441017473
jac: array([54.59815003,  7.3890561 ,  0.36787944, 4.48168907,  0.81873076])
message: 'Optimization terminated successfully.'
nfev: 21
nit: 3
njev: 3
status: 0
success: True
x: array([5.11205896e-11, 1.43702148e-11, 1.00000000e+00, 1.29011784e-11, 0.00000000e+00])


Maybe I'm just missing something, but this optimisation is obviously successful it has weighted the integral most translated to the right.

• Yes in your example the optimization works well. But, when I am using realistic numbers (see edit) I am not able to find a solution in a reasonable timeframe. However, Bondarenko reports in his paper that finding the solution is very fast. Most likely, I have made a mistake in the optimization problem... – William Burknecht Feb 24 at 21:32