4.2 Volatility smile put IV calculations not correct?

Open marcovth opened this issue 1 year ago • 1 comments

Hello ...

Is it possible that something is wrong with the put code of the 4.2 Volatility smile, or is there something I don't understand about the puts?

Thanks for taking a look ...


from FMNM.Parameters import Option_param
from FMNM.Processes import Diffusion_process, Merton_process, VG_process, Heston_process
from FMNM.BS_pricer import BS_pricer
from FMNM.Merton_pricer import Merton_pricer
from FMNM.VG_pricer import VG_pricer
from FMNM.Heston_pricer import Heston_pricer

import numpy as np
import pandas as pd
import scipy as scp
import scipy.stats as ss
import matplotlib.pyplot as plt
import scipy.optimize as scpo
from functools import partial
from itertools import compress
import os
import warnings

warnings.filterwarnings("ignore")

test1=False
if test1:
    spot=2847.60
    otype="put"
    T = 36/365
    S0 = spot
    K = 2835
    strikes = np.arange(800, 4000, 50)  # strike grid

test2=True
if test2:
    spot=100
    otype="put"
    T = 365/365
    S0 = spot
    K = 100
    strikes = np.arange(50, 151, 5)  # strike grid



opt_param_call = Option_param(S0=S0, K=K, T=T, v0=0.04, exercise="European", payoff="call")
opt_param_put = Option_param(S0=S0, K=K, T=T, v0=0.04, exercise="European", payoff="put")

diff_param = Diffusion_process(r=0.1, sig=0.2)
Merton_param = Merton_process(r=0.1, sig=0.1, lam=0.8, muJ=-0.04, sigJ=0.2)
VG_param = VG_process(r=0.1, theta=-0.09, sigma=0.19, kappa=0.6)
Heston_param = Heston_process(mu=0.1, rho=-0.3, sigma=0.6, theta=0.04, kappa=5)

BS_call = BS_pricer(opt_param_call, diff_param)
VG_call = VG_pricer(opt_param_call, VG_param)
Mert_call = Merton_pricer(opt_param_call, Merton_param)
Hest_call = Heston_pricer(opt_param_call, Heston_param)
BS_put = BS_pricer(opt_param_put, diff_param)
VG_put = VG_pricer(opt_param_put, VG_param)
Mert_put = Merton_pricer(opt_param_put, Merton_param)
Hest_put = Heston_pricer(opt_param_put, Heston_param)




BS_prices_call = BS_call.FFT(strikes)
Mert_prices_call = Mert_call.FFT(strikes)
Hest_prices_call = Hest_call.FFT(strikes)
VG_prices_call = VG_call.FFT(strikes)

BS_prices_put = BS_put.FFT(strikes)
Mert_prices_put = Mert_put.FFT(strikes)
Hest_prices_put = Hest_put.FFT(strikes)
VG_prices_put = VG_put.FFT(strikes)

### Closed formula prices
##BS_prices_cl = np.zeros_like(strikes, dtype=float)
##VG_prices_cl = np.zeros_like(strikes, dtype=float)
##Mert_prices_cl = np.zeros_like(strikes, dtype=float)
##Hest_prices_cl = np.zeros_like(strikes, dtype=float)
##
##for i, K in enumerate(strikes):
##    BS.K = K
##    VG.K = K
##    Mert.K = K
##    Hest.K = K
##    BS_prices_cl[i] = BS.closed_formula()
##    VG_prices_cl[i] = VG.Fourier_inversion()
##    Mert_prices_cl[i] = Mert.closed_formula()
##    Hest_prices_cl[i] = Hest.Fourier_inversion()

##print("Closed vs FFT. Total absolute error BS: ", np.linalg.norm(BS_prices - BS_prices_cl, 1))
##print("Closed vs FFT. Total absolute error VG: ", np.linalg.norm(VG_prices - VG_prices_cl, 1))
##print("Closed vs FFT. Total absolute error Merton: ", np.linalg.norm(Mert_prices - Mert_prices_cl, 1))
##print("Closed vs FFT. Total absolute error Heston: ", np.linalg.norm(Hest_prices - Hest_prices_cl, 1))


def implied_volatility(price, S0, K, T, r, payoff="call", method="fsolve", disp=False):
    """Returns Implied volatility
    methods:  fsolve (default) or brent
    """

    def obj_fun(vol):
        if payoff=="call": BS=BS_call
        else: BS=BS_put
        return price - BS.BlackScholes(payoff=payoff, S0=S0, K=K, T=T, r=r, sigma=vol)

    if method == "brent":
        x, r = scpo.brentq(obj_fun, a=1e-15, b=500, full_output=True)
        if r.converged == True:
            return x
    if method == "fsolve":
        X0 = [0.1, 0.5, 1, 3]  # set of initial guess points
        for x0 in X0:
            x, _, solved, _ = scpo.fsolve(obj_fun, x0, full_output=True, xtol=1e-8)
            if solved == 1:
                return x[0]

    if disp == True:
        print("implied_volatility Strike", K)
    return -1


def implied_vol_minimize(price, S0, K, T, r, payoff="call", disp=False):
    """Returns Implied volatility by minimization"""

    n = 2  # must be even

    def obj_fun(vol):
        if payoff=="call": BS=BS_call
        else: BS=BS_put
        return (BS.BlackScholes(payoff=payoff, S0=S0, K=K, T=T, r=r, sigma=vol) - price) ** n

    res = scpo.minimize_scalar(obj_fun, bounds=(1e-15, 8), method="bounded")
    if res.success == True:
        return res.x
    if disp == True:
        print("implied_vol_minimize Strike", K)
    return -1


IV_BS_call = []
IV_VG_call = []
IV_Mert_call = []
IV_Hest_call = []

IV_BS_put = []
IV_VG_put = []
IV_Mert_put = []
IV_Hest_put = []
for i in range(len(strikes)):
    IV_BS_call.append(implied_volatility(BS_prices_call[i], S0=spot, K=strikes[i], T=T, r=0.1))
    IV_VG_call.append(implied_volatility(VG_prices_call[i], S0=spot, K=strikes[i], T=T, r=0.1))
    IV_Mert_call.append(implied_volatility(Mert_prices_call[i], S0=spot, K=strikes[i], T=T, r=0.1))
    IV_Hest_call.append(implied_volatility(Hest_prices_call[i], S0=spot, K=strikes[i], T=T, r=0.1))

    IV_BS_put.append(implied_volatility(BS_prices_put[i], S0=spot, K=strikes[i], T=T, r=0.1))
    IV_VG_put.append(implied_volatility(VG_prices_put[i], S0=spot, K=strikes[i], T=T, r=0.1))
    IV_Mert_put.append(implied_volatility(Mert_prices_put[i], S0=spot, K=strikes[i], T=T, r=0.1))
    IV_Hest_put.append(implied_volatility(Hest_prices_put[i], S0=spot, K=strikes[i], T=T, r=0.1))

##IV_BS_m = []
##IV_VG_m = []
##IV_Mert_m = []
##IV_Hest_m = []
##for i in range(len(strikes)):
##    #IV_BS_m.append(implied_vol_minimize(BS_prices[i], S0=spot, K=strikes[i], T=T, r=0.1))
##    #IV_VG_m.append(implied_vol_minimize(VG_prices[i], S0=spot, K=strikes[i], T=T, r=0.1))
##    IV_Mert_m.append(implied_vol_minimize(Mert_prices[i], S0=spot, K=strikes[i], T=T, r=0.1))
##    #IV_Hest_m.append(implied_vol_minimize(Hest_prices[i], S0=spot, K=strikes[i], T=T, r=0.1))


##print(
##    " Are the IV values obtained by the methods above equal? ",
##    np.allclose(np.array(IV_BS), np.array(IV_BS_m))
##    & np.allclose(np.array(IV_VG), np.array(IV_VG_m))
##    & np.allclose(np.array(IV_Mert), np.array(IV_Mert_m))
##    & np.allclose(np.array(IV_Hest), np.array(IV_Hest_m)),
##)

fig = plt.figure(figsize=(16, 4))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

ax1.plot(strikes, BS_prices_call, label="BS_call")
ax1.plot(strikes, VG_prices_call, label="VG_call")
ax1.plot(strikes, Mert_prices_call, label="Mert_call")
ax1.plot(strikes, Hest_prices_call, label="Hest_call")

ax1.plot(strikes, BS_prices_put, label="BS_put")
ax1.plot(strikes, VG_prices_put, label="VG_put")
ax1.plot(strikes, Mert_prices_put, label="Mert_put")
ax1.plot(strikes, Hest_prices_put, label="Hest_put")


#ax1.set_xlim([80, 150])
#ax1.set_ylim([0, 30])
ax1.set_title("Comparison of prices")
ax1.set_xlabel("Strike")
ax1.set_ylabel("Price")

ax2.plot(strikes, IV_BS_call, label="IV_BS_call")
ax2.plot(strikes, IV_VG_call, label="IV_VG_call")
ax2.plot(strikes, IV_Mert_call, label="IV_Mert_call")
ax2.plot(strikes, IV_Hest_call, label="IV_Hest_call")

ax2.plot(strikes, IV_BS_put, label="IV_BS_put")
ax2.plot(strikes, IV_VG_put, label="IV_VG_put")
ax2.plot(strikes, IV_Mert_put, label="IV_Mert_put")
ax2.plot(strikes, IV_Hest_put, label="IV_Hest_put")



ax2.set_title("Comparison of Implied volatilities")
ax2.set_xlabel("Strike")
ax2.set_ylabel("Imp Vol")
ax1.legend()
ax2.legend()
plt.show()

Apr 07 '25 21:04 marcovth

Hi @marcovth Can you be more precise about the error you found? There is no need to copy/paste irrelevant code. From a first look, I can see you have convergence problems with the implied volatility function. Consider that:

when the implied volatility solver fails, it returns -1.
The integrals in the FFT method may not always converge because they were set for the cases used in the notebook.

Apr 08 '25 08:04 cantaro86