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