Implied Volatility: Newton-Raphson and Bisection Method

Newton-Raphson Method

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

def bsm_option_price(S, K, T, r, sigma, option_type='call'):
    d1 = (math.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * math.sqrt(T))
    d2 = d1 - sigma * math.sqrt(T)
    if option_type == 'call':
        return S * norm.cdf(d1) - K * math.exp(-r * T) * norm.cdf(d2)
    elif option_type == 'put':
        return K * math.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)

def calculate_implied_volatility(option_price, S, K, T, r, option_type='call'):
    epsilon = 1e-6
    sigma = 0.5  # Initial guess for volatility
    max_iterations = 100
    iterations = 0

    # Lists to store iteration results for plotting
    iteration_values = []
    iv_values = []

    while True:
        d1 = (math.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * math.sqrt(T))
        vega = S * norm.pdf(d1) * math.sqrt(T)
        option_price_estimate = bsm_option_price(S, K, T, r, sigma, option_type)

        iteration_values.append(iterations)
        iv_values.append(sigma)

        error = option_price - option_price_estimate
        if abs(error) < epsilon or iterations >= max_iterations:
            break

        sigma += error / vega
        iterations += 1

    return sigma, iteration_values, iv_values

# Example values for option parameters
S = 100  # Underlying price
K = 100  # Strike price
T = 1  # Time to expiration
r = 0.02  # Risk-free rate
option_price = 10  # Market price of the option
option_type = 'call'  # 'call' or 'put'

implied_volatility, iteration_values, iv_values = calculate_implied_volatility(option_price, S, K, T, r, option_type)

print(f"Implied Volatility: {implied_volatility:.6f}")

# Display iteration results
print("Iterations for Implied Volatility Calculation:")
print("Iteration |  Implied Volatility  |  Error")
print("-----------------------------------------")
for i, sigma in enumerate(iv_values):
    error = option_price - bsm_option_price(S, K, T, r, sigma, option_type)
    print(f"   {i+1:2}     |       {sigma:.6f}       |  {error:.8f}")

# Plotting the iteration of IV values
plt.figure(figsize=(10, 6))
plt.plot(iteration_values, iv_values, marker='o', linestyle='-')
plt.xlabel('Iterations')
plt.ylabel('Implied Volatility')
plt.title('Newton-Raphson Method: Implied Volatility Iteration')
plt.grid(True)
plt.show()

Output

---

Bisection Method

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

def bsm_option_price(S, K, T, r, sigma, option_type='call'):
    d1 = (math.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * math.sqrt(T))
    d2 = d1 - sigma * math.sqrt(T)
    if option_type == 'call':
        return S * norm.cdf(d1) - K * math.exp(-r * T) * norm.cdf(d2)
    elif option_type == 'put':
        return K * math.exp(-r * T) * norm.cdf(-d2) - S * norm.cdf(-d1)

def calculate_implied_volatility_bisection(option_price, S, K, T, r, option_type='call'):
    epsilon = 1e-6
    low_volatility = 0.01
    high_volatility = 3.0
    max_iterations = 100
    iterations = 0

    # Lists to store iteration results for plotting
    iteration_values = []
    iv_values = []

    while iterations < max_iterations:
        mid_volatility = (low_volatility + high_volatility) / 2.0
        option_price_estimate = bsm_option_price(S, K, T, r, mid_volatility, option_type)

        iteration_values.append(iterations)
        iv_values.append(mid_volatility)

        error = option_price - option_price_estimate

        if abs(error) < epsilon:
            break
        elif error < 0:
            high_volatility = mid_volatility
        else:
            low_volatility = mid_volatility

        iterations += 1

    return mid_volatility, iteration_values, iv_values

# Example values for option parameters
S = 100  # Underlying price
K = 100  # Strike price
T = 1  # Time to expiration
r = 0.02  # Risk-free rate
option_price = 10  # Market price of the option
option_type = 'call'  # 'call' or 'put'

implied_volatility, iteration_values, iv_values = calculate_implied_volatility_bisection(option_price, S, K, T, r, option_type)

print(f"Implied Volatility (Bisection Method): {implied_volatility:.6f}")

# Display iteration results in a table
print("Iterations for Implied Volatility Calculation (Bisection Method):")
print("Iteration |  Implied Volatility  |  Option Price Error")
print("----------------------------------------------------------")
for i, sigma in enumerate(iv_values):
    error = option_price - bsm_option_price(S, K, T, r, sigma, option_type)
    print(f"   {i+1:2}     |       {sigma:.6f}       |  {error:.8f}")

# Plotting the iteration of IV values
plt.figure(figsize=(10, 6))
plt.plot(iteration_values, iv_values, marker='o', linestyle='-')
plt.xlabel('Iterations')
plt.ylabel('Implied Volatility')
plt.title('Bisection Method: Implied Volatility Iteration')
plt.grid(True)
plt.show()

Output