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

Question

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.integrate import dblquad
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
    return total

# 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))
    return total

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

Answer 1

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}} \;. $$

Then your optimisation problem becomes:

$$\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...

Answer 2

Having read your comments I'm not sure where the implementation issue is:

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
from scipy.integrate import dblquad

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.