My question is closely related to the answer of @LocalVolatility and his blogpost. I am trying to reproduce his first figure and I am struggling with the implied volatility. With the help of

$$ f(S) = \sum^n_{i=1}p_i\frac{1}{\sqrt{2\pi}\sigma_i}exp\{-\frac{(S-\mu_i)^2}{2\sigma_i^2}\} $$

I know the probability density function. I also know that the price of a call is equal to

$$ C=e^{-rT}\int_0^\infty(S-K)^+p(S)dS $$

I write the following Python code

def call(r, T, S, K, f):
    inner = np.zeros_like(S)
    itm = S > K
    inner[itm] = S[itm] - K

    return np.exp(-r * T) * np.sum(inner * f)

For $K = \{90, 91, ..., 110\}$ I calculate the call price and use this price to infer the IV. My plot looks as follows:

enter image description here

The shape of my implied volatility smile is different. What am I doing wrong? Intuitively, I would say it makes sense that the IV of the C110 is zero as we don't expect the underlying to trade above 110.

Full Python code after tips of @LocalVolatility

import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt

def f(mu, sigma, x):
    return np.exp((-1 * (x - mu) ** 2) / (2 * sigma ** 2)) / (np.sqrt(2 * np.pi) * sigma)

def call(r, T, S, K, pdf, dx):
    inner = np.zeros_like(S)
    itm = S > K
    inner[itm] = S[itm] - K

    return np.exp(-r * T) * np.sum(inner * pdf) * dx

def call_bs(S, K, T, r, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

    return S * ss.norm.cdf(d1, 0, 1) - K * np.exp(-r * T) * ss.norm.cdf(d2, 0, 1)

def vega_bs(S, K, T, r, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

    return S * ss.norm.pdf(d1) * np.sqrt(T)

def iv(c, S, K, T, r):
    N_ITERATIONS = 200
    THRESH = 1.0E-5
    sigma = 0.25

    for _ in range(N_ITERATIONS):
        p = call_bs(S, K, T, r, sigma)
        v = vega_bs(S, K, T, r, sigma)
        diff = c - p

        if abs(diff) < THRESH:
            return sigma

        sigma = sigma + diff / v

    return sigma

# input
S = 100
dx = 0.1
n = 2
p_up = 0.5
p_down = 1 - p_up
mu_up = 0.05
sigma_up = 0.02
mu_down = -0.05
sigma_down = 0.02
T = 7 / 255
r = 0

if __name__ == "__main__":
    # make integration zone wider
    x = np.arange(S - 25, S + 25, dx)

    # use log returns; not prices
    f_up = p_up * f(mu=mu_up, sigma=sigma_up, x=np.log(x / S))
    f_down = p_down * f(mu=mu_down, sigma=sigma_down, x=np.log(x / S))
    implied_density = (f_up + f_down) / S
    area_under_implied_density = np.trapz(implied_density, dx=dx) # 1.0014505305168542

    ax1 = plt.subplot()
    ax1.plot(x, implied_density, c='blue', label='implied density')
    ax1.set_ylabel('implied density')
    ax1.legend(loc='upper left')
    ax1.set_ylim([0, 0.12])
    ax2 = ax1.twinx()

    ivs = []

    for _ in x:
        c = call(r=r, T=T, S=x, K=_, pdf=implied_density, dx=dx)
        implied_vol = iv(c=c, S=S, K=_, T=T, r=r)

    ax2.plot(x, ivs, c='green', label='implied volatility')
    ax2.set_ylabel('implied volatility')

    plt.grid(axis='both', linestyle='--')
    plt.xlim([90, 110])
    ax2.legend(loc='upper right')

The plot now looks as follows:

enter image description here

I am still doing something wrong for the ITM-options. Any suggestions?

  • 2
    $\begingroup$ Can you share what you're using to get the implied volatility? $\endgroup$
    – KT8
    Commented Feb 19 at 13:11
  • 1
    $\begingroup$ (1) Can you give your prices for a few options please. E.g. strike 90, 95, 100, 105, 110? (2) What is the grid for S that you use in pricing? $\endgroup$ Commented Feb 26 at 14:15
  • $\begingroup$ 90: 99.4219 95: 53.5350 100: 24.7310 105: 3.8241 110: 0 This is when I am dividing S up in increments of 0.1 $\endgroup$
    – HJA24
    Commented Feb 26 at 14:34

2 Answers 2


This is not an answer to the question but too long for a comment. There are a bunch of issues with your code. In roughly decreasing order of importance:

  1. Sanity-checking your prices should immediately tell you that something is not right there. Your strike 90 USD call is worth 99.42 USD?

  2. The main mistake in your pricing function is that you forget to normalize by the grid size dx when integrating. I.e. replace np.sum(inner * pdf) by np.sum(inner * pdf) * dx. Now things already look a bit better.

  3. You are imposing a normal distribution on prices while you should model the log-returns as normally distributed. When I talk about a normally distributed jump size in my post, then this is always in log-return space. I.e. you need to replace your density to be a log-normal density.

  4. When you model the log-returns as normal, you will need to do a martingale correction to ensure the forward doesn't change around the jump. I.e. you need to ensure that $p e^{\mu_+} + (1 - p) e^{\mu_-} = 1$ by shifting one of the three accordingly. You cannot impose all of them!

  5. Generally, you should choose your integration range to be wider and de-couple your spot and strike grids. You e.g. implicitly attribute zero probability to terminal spots below 90 USD while it is only 2.5 standard deviations down from the approximate mean in one of the two states. The integral of your density over the spot grid is "only" 0.9942, i.e. you're loosing some mass.

  6. Your code will be broken when you actually want to use non-zero rates because you don't account for this in your density / your terminal forward will be wrong.

As a general recommendation: You need to sanity-check your own work better. Don't jump all the way to the result but make sure some basic relationships hold.

  • E.g. one thing I always test is that the price of a zero-strike call option is equal to the discounted forward.

  • You could also check that the the price with both jump means being zero is equal to the Black-Scholes price when matching the standard deviation.

These tests would have already failed in your code, drawing your attention to the pricing part.

  • $\begingroup$ Thank you for your tips! Is step 4 really necessary because when I calculate the area of the implied density it is +-1? In addition, do you see why the implied skew for the ITM-options is different? $\endgroup$
    – HJA24
    Commented Mar 5 at 13:47
  • 1
    $\begingroup$ Any properly constructed density integrates to one but that doesn't mean the corresponding expectation is correct. You know by no arbitrage that your expectation has to be equal to the forward price both in the left and right limit at the jump time. $\endgroup$ Commented Mar 5 at 13:49
  • 1
    $\begingroup$ Your skew being wrong on the upside means something else is wrong still in your pricing. I didn't have an in-depth look but my best guess is your forward is still wrong. Did you do the sanity-checks that I suggested? Probably missing the $-\sigma^2 / 2$ correction term in the log-normal mean. $\endgroup$ Commented Mar 5 at 14:04
  • $\begingroup$ So, the Newton IV's method is superfluous? $\endgroup$
    – HJA24
    Commented Mar 7 at 9:12
  • $\begingroup$ In your blogpost you mention: "The base dynamics follow a constant coefficient geometric Brownian motion with $\sigma = 15% $" Where does this come back in the methodology? $\endgroup$
    – HJA24
    Commented Mar 7 at 9:19

It looks like your approach to calculating IV is generally on the right track. However, to ensure accuracy and address any discrepancies in the implied volatility smile you're observing, I would refine methodology in two areas:

1. Integration for Call Price Calculation: The method you've described for calculating the call price using a probability density function (PDF) (f(S)) seems to manually sum over discrete values. For a more accurate calculation, especially if (S) is not densely sampled, I would suggest integrating over a continuous range of (S) using libraries like scipy:

import numpy as np
import scipy.stats as ss
import scipy.integrate as integrate

def f(S, pis, mus, sigmas):
    """Probability density function for a mixture of normals."""
    result = 0
    for pi, mu, sigma in zip(pis, mus, sigmas):
        result += pi / (np.sqrt(2 * np.pi) * sigma) * np.exp(-((S - mu)**2) / (2 * sigma**2))
    return result

def call_price(r, T, K, pis, mus, sigmas):
    """Calculate call price using numerical integration."""
    f_integrated = lambda S: np.exp(-r * T) * (S - K) * f(S, pis, mus, sigmas)
    price, _ = integrate.quad(f_integrated, K, np.inf)
    return price

# Example parameters for f(S)
pis = [0.5, 0.5]  # Weights
mus = [-0.05, 0.05]  # Means
sigmas = [0.02, 0.02]  # Standard deviations

# Call price calculation example
r = 0.00  # Risk-free rate
T = 1/52  # Time to maturity in years
K = 100  # Strike price

# Adjust pis, mus, sigmas according to your specifics
price_example = call_price(r, T, K, pis, mus, sigmas)
print(f"Example call price: {price_example}")

2. Implied Volatility Calculation: Your iterative approach to finding the IV using the Newton-Raphson method seems appropriate. The initial guess and the iteration limit seem reasonable to me, but I’d refine the calculation of the option price (price) and vega (vega) are accurately calculating the Black-Scholes price and vega, respectively (see Haug and Hull for detailed explanations). The key is to ensure that the call price calculation (call_price) matches the market or theoretical prices you're using to infer IV:

import numpy as np
import scipy.stats as ss

def black_scholes_price(S, K, T, r, sigma):
    """Calculate Black-Scholes option price."""
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    call_price = S * ss.norm.cdf(d1, 0.0, 1.0) - K * np.exp(-r * T) * ss.norm.cdf(d2, 0.0, 1.0)
    return call_price

def black_scholes_vega(S, K, T, r, sigma):
    """Calculate Black-Scholes Vega, the derivative of the price with respect to volatility."""
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    vega = S * ss.norm.pdf(d1) * np.sqrt(T)
    return vega

def find_implied_volatility(target_price, S, K, T, r):
    """Find the implied volatility given a target option price."""
    PRECISION = 1.0e-5
    sigma_guess = 0.5
    for _ in range(MAX_ITERATIONS):
        price = black_scholes_price(S, K, T, r, sigma_guess)
        vega = black_scholes_vega(S, K, T, r, sigma_guess)
        price_diff = target_price - price

        if abs(price_diff) < PRECISION:
            return sigma_guess

        sigma_guess += price_diff / vega

        # Prevent negative volatility
        if sigma_guess < 0:
            sigma_guess = PRECISION

    return sigma_guess

# Example usage:
S = 100  # Underlying asset price
K = 100  # Strike price
T = 1/52  # Time to expiration in years (1 week)
r = 0.00  # Risk-free interest rate
target_price = 2.5  # Target call price from market or model

implied_volatility = find_implied_volatility(target_price, S, K, T, r)
print(f"Implied Volatility: {implied_volatility}")


Haug, E. G. (2007). The Complete Guide to Option Pricing Formulas (2nd ed.). McGraw-Hill Education.

Hull, J. C. (2017). Options, Futures, and Other Derivatives (10th ed.). Pearson.

  • $\begingroup$ Thanks for your time, but your suggestions do not lead to a significant change of the volatility smile. You're basically repeating what I already stated that I was doing... So, my question remains unanswered. $\endgroup$
    – HJA24
    Commented Feb 20 at 8:11
  • $\begingroup$ Due to the high dynamic range of option prices, I tend to prefer halting conditions on implied volatility solvers that test the change in volatility parameter, rather than the price difference. Nice code (upvoted), though, @Theo $\endgroup$
    – Brian B
    Commented Feb 26 at 19:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.