# Lognormal-mixture dynamics and calibration to market volatility smiles

Can someone assist me in replicating the code and results from page 11, Figure 3 of the paper 'Lognormal-mixture dynamics and calibration to market volatility smiles' by Damiano Brigo, Fabio Mercurio, and Francesco Rapisarda, available at https://www.ma.imperial.ac.uk/~dbrigo/lognsmile.pdf? I have posted my code attempt, but it does not perfectly match the one in the paper. Any help or guidance on achieving this would be greatly appreciated.

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

class LognormalMixture(object):
"""
*
* Title: Lognormal-mixture dynamics and calibration
* to market volatility smiles.
* From: Damiano Brigo, Fabio Mercurio, Francesco Rapisarda.
*
"""
class Shifted(object):
"""
*
* From chapter 4: Shifting the overall distribution.
* Quote: "We now show how to construct an even more general model
* by shifting the process (11), while preserving the correct drift.
* Precisely, we assume that the new asset-price process A0 is obtained
* through the following affine transformation of the process S ...".
*
"""
def _sigma_0 \
(
lambda_: np.ndarray, eta_: np.ndarray, t: float, alpha: float
) -> float:
"""
*
* \sigma(0) is the ATM-forward implied volatility.
*
* @param array $$lambda_ - weights of each lognormal dist. * @param array$$eta_    - volatilities standardized by time
* @param float $$t - maturity * @param float$$alpha   -
*
"""
return 2/np.sqrt(t) * norm.ppf((1 - alpha) * np.sum(lambda_ * \
norm.cdf(0.5 * eta_ * np.sqrt(t))) + 0.5 * alpha)

def _sigma_m \
(
k: np.ndarray, s0: float, mu: float, lambda_: np.ndarray,
eta_: np.ndarray, t: float, alpha: float
) -> float:
"""
*
* Shifting the overall distribution:
* Black-Scholes volatility derived from array of log-moneyness.
*
* @param array $$m - log moneyness * @param array$$lambda_ - weights of each lognormal dist.
* @param array $$eta_ - volatilities standardized by time * @param float$$t       - maturity
* @param float \$alpha   -
*
"""
# calculte log-moneyness
k_a = k - s0 * alpha * np.exp(mu*t)
a0 = s0 * (1 - alpha)
m = np.log(a0/k_a) + mu * t

# atm-forward implied volatilities
sigma_0: float = LognormalMixture.Shifted._sigma_0(
lambda_, eta_, t, alpha)

# black-scholes implied-volatilities
sigma_m = sigma_0
sigma_m += alpha * m * (np.sum(lambda_ * \
norm.cdf(-0.5 * eta_ * np.sqrt(t))) - 0.5) * (np.sqrt(t) / \
np.sqrt(np.pi*2) * np.exp(-0.125 * sigma_0**2 * t))**-1
return sigma_m + 0.5 * ((t*(1-alpha))**-1 * np.sum(lambda_/eta_ * \
np.exp(0.125*t*(sigma_0**2-eta_**2))) - (sigma_0*t)**-1 + \
0.25 * alpha**2 * sigma_0 * t * ((np.sum(lambda_ * norm.cdf(
-0.5 * eta_ * np.sqrt(t))) - 0.5)/ (np.sqrt(t)/np.sqrt(
np.pi*2) * np.exp(-0.125 * sigma_0**2 * t)))**2) * m**2

lnm = LognormalMixture.Shifted
k_ = np.arange(80,121,0.5)
lambda_ = np.array([0.6,0.4])
s0 = 100
t = 2
plt.figure(figsize=(11,7))

# sp1
plt.subplot(121)
plt.grid()
plt.plot(k_, lnm._sigma_m(k=k_, s0=s0, mu=0.05, lambda_=lambda_,
eta_=np.array([0.35,0.1]),  t=t, alpha=0.0), label="α=0.0", c="b")
plt.plot(k_, lnm._sigma_m(k=k_, s0=s0, mu=0.05, lambda_=lambda_,
eta_=np.array([0.35,0.1]),  t=t, alpha=-0.2), label="α=-0.2", c="g", ls="--")
plt.plot(k_, lnm._sigma_m(k=k_, s0=s0, mu=0.05, lambda_=lambda_,
eta_=np.array([0.35,0.1]),  t=t, alpha=-0.4), label="α=-0.4", c="r", ls=":")
plt.legend()

# sp2
plt.subplot(122)
plt.grid()
plt.plot(k_, lnm._sigma_m(k=k_, s0=s0, mu=0.05, lambda_=lambda_,
eta_=np.array([0.35,0.1]),  t=t, alpha=0.0), label="α=0.0", c="b")
plt.plot(k_, lnm._sigma_m(k=k_, s0=s0, mu=0.05, lambda_=lambda_,
eta_=np.array([0.1099,0.3553]),  t=t, alpha=-0.2), label="α=-0.2", c="g", ls="--")
plt.plot(k_, lnm._sigma_m(k=k_, s0=s0, mu=0.05, lambda_=lambda_,
eta_=np.array([0.09809,0.2979]),  t=t, alpha=-0.4), label="α=-0.4", c="r", ls=":")
plt.legend()