from rnd_cpp import *
from scipy.optimize import minimize
import numpy as np
import time
import warnings
warnings.simplefilter("ignore")
def ExtractMlnDensity(r,
y,
te,
s0,
market_calls,
call_strikes,
market_puts,
put_strikes,
call_weights = 1,
put_weights = 1,
lambda_ = 1,
hessian_flag = False,
initial_values = np.repeat(np.NaN, 5),
cl = {'maxit': 10000}):
if np.sum(np.isnan(initial_values)) >= 1:
initial_values = extract_mln_density_grid(r,
y,
te,
s0,
market_calls,
call_strikes,
market_puts,
put_strikes)
optim_obj = minimize(fun = MlnObjective,
x0 = initial_values,
method = 'Nelder-Mead',
bounds = [(0, 1), (0, 1e10), (0, 1e10), (0, 1e10), (0, 1e10)],
args = (r, y, te, s0, market_calls, call_strikes, market_puts, put_strikes),
options={'ftol': 1e-04, 'maxiter': 10000, 'disp': False})
alpha_1 = optim_obj.x[0]
meanlog_1 = optim_obj.x[1]
meanlog_2 = optim_obj.x[2]
sdlog_1 = optim_obj.x[3]
sdlog_2 = optim_obj.x[4]
converge_result = optim_obj.status
out = {'alpha_1': alpha_1, 'meanlog_1': meanlog_1, 'meanlog_2': meanlog_2, 'sdlog_1': sdlog_1, 'sdlog_2': sdlog_2}
return out
def MlnObjective(theta,
r,
y,
te,
s0,
market_calls,
call_strikes,
market_puts,
put_strikes,
call_weights = 1,
put_weights = 1,
lambda_ = 1):
return mln_objective(theta,
r,
y,
te,
s0,
market_calls,
call_strikes,
market_puts,
put_strikes,
call_weights,
put_weights,
lambda_)
def PriceMlnOption(r,
te,
y,
k,
alpha_1,
meanlog_1,
meanlog_2,
sdlog_1,
sdlog_2):
return price_mln_option(r,
te,
y,
k,
alpha_1,
meanlog_1,
meanlog_2,
sdlog_1,
sdlog_2)
# def dmln(x, alpha_1, meanlog_1, meanlog_2, sdlog_1, sdlog_2):
# return dmln(x,
# alpha_1,
# meanlog_1,
# meanlog_2,
# sdlog_1,
# sdlog_2)
if __name__ == '__main__':
import matplotlib.pyplot as plt
#
# Create some calls and puts based on mln and
# see if we can extract the correct values.
#
r = 0.05
y = 0.02
te = 60/365
meanlog_1 = 6.8
meanlog_2 = 6.95
sdlog_1 = 0.065
sdlog_2 = 0.055
alpha_1 = 0.4
call_strikes = np.arange(800, 1210, 10)
market_calls = PriceMlnOption(r = r, y = y, te = te, k = call_strikes, alpha_1 = alpha_1, meanlog_1 = meanlog_1, meanlog_2 = meanlog_2, sdlog_1 = sdlog_1, sdlog_2 = sdlog_2)['call']
s0 = PriceMlnOption(r = r, y = y, te = te, k = call_strikes, alpha_1 = alpha_1, meanlog_1 = meanlog_1, meanlog_2 = meanlog_2, sdlog_1 = sdlog_1, sdlog_2 = sdlog_2)['s0']
put_strikes = np.arange(800, 1210, 10)
market_puts = PriceMlnOption(r = r, y = y, te = te, k = call_strikes, alpha_1 = alpha_1, meanlog_1 = meanlog_1, meanlog_2 = meanlog_2, sdlog_1 = sdlog_1, sdlog_2 = sdlog_2)['put']
t1 = time.time()
print(ExtractMlnDensity(r,
y,
te,
s0[0],
market_calls,
call_strikes.astype(float).tolist(),
market_puts,
put_strikes.astype(float).tolist()))
print(time.time() - t1)
#
# The mln objective function should be close to zero.
# The weights are automatically set to 1.
#
r = 0.05
te = 60/365
y = 0.02
meanlog_1 = 6.8
meanlog_2 = 6.95
sdlog_1 = 0.065
sdlog_2 = 0.055
alpha_1 = 0.4
# This is the current price implied by parameter values:
s0 = 981.8815
call_strikes = np.arange(800, 1210, 10)
market_calls = PriceMlnOption(r=r, y = y, te = te, k = call_strikes,
alpha_1 = alpha_1, meanlog_1 = meanlog_1, meanlog_2 = meanlog_2,
sdlog_1 = sdlog_1, sdlog_2 = sdlog_2)['call']
put_strikes = np.arange(800, 1210, 10)
market_puts = PriceMlnOption(r=r, y = y, te = te, k = call_strikes,
alpha_1 = alpha_1, meanlog_1 = meanlog_1, meanlog_2 = meanlog_2,
sdlog_1 = sdlog_1, sdlog_2 = sdlog_2)['put']
MlnObjective(theta=[alpha_1,meanlog_1, meanlog_2 , sdlog_1, sdlog_2],
r = r, y = y, te = te, s0 = s0,
market_calls = market_calls, call_strikes = call_strikes,
market_puts = market_puts, put_strikes = put_strikes, lambda_ = 1)
#
# Try out a range of options
#
r = 0.05
te = 60/365
k = np.arange(700, 1301, 1)
y = 0.02
meanlog_1 = 6.80
meanlog_2 = 6.95
sdlog_1 = 0.065
sdlog_2 = 0.055
alpha_1 = 0.4
mln_prices = PriceMlnOption(r = r, y = y, te = te, k = k, alpha_1 = alpha_1,
meanlog_1 = meanlog_1, meanlog_2 = meanlog_2, sdlog_1 = sdlog_1, sdlog_2 = sdlog_2)
plt.scatter(x = k, y = mln_prices['call'], c = 'black')
plt.scatter(x = k, y = mln_prices['put'], c = 'red')
plt.grid()
plt.show()
#
# A bimodal risk neutral density!
#
mln_alpha_1 = 0.4
mln_meanlog_1 = 6.3
mln_meanlog_2 = 6.5
mln_sdlog_1 = 0.08
mln_sdlog_2 = 0.06
k = np.arange(300, 900)
dx = dmln(x = k, alpha_1 = mln_alpha_1, meanlog_1 = mln_meanlog_1,
meanlog_2 = mln_meanlog_2,
sdlog_1 = mln_sdlog_1, sdlog_2 = mln_sdlog_2)
plt.plot(k, dx)
plt.grid()
