Return to Question

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 to 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(500, 800, num=n)

rnd = np.empty((1, n))
for d in range(0, j):
    rnd += norm.pdf(s,z[d],h) * a_sol[d]
    
plt.figure()
plt.plot(s, rnd[0])

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)

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.

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 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 to 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(500, 800, num=n)

rnd = np.empty((1, n))
for d in range(0, j):
    rnd += norm.pdf(s,z[d],h) * a_sol[d]
    
plt.figure()
plt.plot(s, rnd[0])

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.