使用正卷積近似估計風險中性密度 - Python
我試圖通過正卷積近似來估計風險中性密度（由 Bondarenko 2002 介紹：https ://papers.ssrn.com/sol3/papers.cfm?abstract_id=375781 ）。我目前正在努力在 Python 中實現計算算法來解決以下優化問題：
$$ \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 $$
在哪裡： $ P_i $ 是觀察到的有行使價的看跌期權價格 $ x_i $ , $ \phi $ 參考重新調整的標準正態分佈和 $ z_j $ 表示最小和最大觀察到的罷工之間等距網格上的點 $ x $ .
很簡單，應該使用標準的二次程序來解決問題。但是，我不知道如何處理 $ a $ 是在一個函式內 $ u $ ，它本身在一個二重積分內。
是否有人已經實現了正卷積近似來估計 Python 中的風險中性密度？
否則，有人可以告訴我如何用雙積分編寫優化問題，例如：
$$ \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 $$
謝謝您的幫助！
編輯
*更新：*感謝 Attack68 的評論和回答。我能夠實現以下程式碼：
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)它返回以下警告並且難以收斂：
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)在堆棧溢出文章之後，我將嘗試使用 nquad 而不是 dblquad。我將發布進一步的進展。
編輯 2 更新：使用 Attack68 的第二個答案的見解，我能夠以“有效”的方式估計 RND（可能可以進一步改進）：
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)
再看一遍，你的實現可能不是很低效，因為； $$ \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);, $$
在哪裡 $ \Phi(y) $ 是轉換後的累積正態分佈函式 $ y $ . 此功能已經存在
scipy.stats.norm.cdf並已優化。那麼你有; $$ \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 ;. $$
這些 $ K_j $ 只是從的角度來看的常數 $ a_j $ 優化所以可以預先計算，而且既然你有
scipy.stats.norm.cdf它，你似乎只需要使用quad而不是dblquad. 您可能還想知道對於標準的普通 cdf（請參閱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 $$
通過更多的分析，我想知道是否不能進一步簡化……