Pricing Options from an Implied Distribution and Recovering Volatility Skew

# Infer implied volatility skew/smile from implied distribution


# Infer implied volatility skew/smile from implied distribution












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:

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)
        ivs.append(implied_vol)

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

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

The plot now looks as follows:

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

## Answer by LocalVolatility (score 2, accepted)

https://quant.stackexchange.com/a/78571

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:

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

- 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.

- 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.

- 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!

- 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.

- 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.

## Answer by Theo (score 1)

https://quant.stackexchange.com/a/78397

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."""
    MAX_ITERATIONS = 200
    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}")
```

References

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.