Estimación de Riesgo-Neutro Utilizando Densidades Positivo de Convolución de Aproximación - Python

Question

Estoy tratando de estimar el riesgo de densidad neutra positivo a través de la convolución de la aproximación (introducido por Bondarenko 2002: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=375781). Actualmente estoy luchando para implementar un algoritmo computacional en Python solucionar el siguiente problema de optimización:

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

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

donde: $P_i$ es el observado ponga el precio con la huelga de $x_i$, $\phi$ referirse a un reescalado distribución normal estándar y $z_j$ representa los puntos en un igualmente espaciado de cuadrícula entre la más pequeña y la más grande observado huelgas de $x$.

Es fácil que uno debe usar un estándar cuadrática programa para resolver el problema. Sin embargo, no sé cómo manejar el hecho de que el $a$'s están dentro de una función de $u$, que a su vez está dentro de una integral doble.

¿Alguien ya implementado positivo de convolución aproximación para estimar el riesgo de densidad neutra en Python?

De lo contrario, es posible que alguien me muestre cómo el código de un problema de optimización con una integral doble, como por ejemplo:

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

Gracias por la ayuda!

---

EDITAR

Actualización: Gracias a los comentarios y la respuesta de Attack68. Yo era capaz de implementar el siguiente código:

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)

que devuelve la siguiente advertencia y tiene dificultad para convergen:

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)

Tras un desbordamiento de pila post voy a tratar de utilizar nquad en lugar de dblquad. Voy a publicar más progreso.

---

EDIT 2 Actualización: el Uso de las pistas de Attack68 la segunda respuesta, que fue capaz de estimar la RND en una "eficaz" de forma (probablemente puede ser mejorado):

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

Después de haber leído sus comentarios no estoy seguro de que la cuestión de la ejecución es:

scipy.integrate.dblquad resuelve la siguiente integral iterada problema:

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

Si me re-expresar sus variables para ser coherente con la scipy documentación, y supongo que he cambiado la integral finita:

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

Entonces no sería tan simple como:

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]

Y para involucrar a la optimización usted puede intentar:

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

A lo mejor me estoy perdiendo algo, pero esta optimización es, obviamente, el éxito se ha ponderado la integral más traducido a la derecha.

Answer 2

Mirando un segundo momento, su aplicación no es probablemente muy ineficiente, ya que para; $$ \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)\;, $$

donde $\Phi(x)$ es la transformada acumulativa normal de la función de distribución de $y$. Esta función ya existe como scipy.stats.norm.cdf y es la óptima.

Entonces usted tiene; $$ \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 \;.$$

Estos $K_j$ sólo son constantes desde el punto de vista de la $a_j$ optimización por lo que puede ser precalculadas, y ya que usted scipy.stats.norm.cdf parece que usted sólo tiene que utilizar quad e no dblquad. Usted también puede estar interesado en saber que para que una normal estándar cdf (ver 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}} \;. $$

Luego de su optimización problema se convierte en:

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

Con un poco más de análisis me pregunto si eso no puede ser simplificado aún más...