Exploring Greeks – Gamma (Γ) and Theta (Θ) with Python - Simplify Quant

Relationship of Gamma (Γ) with different expiries

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

# Function to calculate call gamma
def calculate_call_gamma(S, K, r, T, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    gamma = np.exp(-r * T) * norm.pdf(d1) / (S * sigma * np.sqrt(T))
    return gamma

# Parameters
S = 100  # Current stock price
K_range = np.arange(80, 121, 5)  # Strike price range at intervals of 5
sigma = 0.3  # Implied volatility
r = 0.05  # Risk-free rate
expiry_periods = [1/12, 3/12, 6/12]  # Expiry periods in years (1 month, 3 months, 6 months)

# Calculate call gamma for different strike prices and expiry periods
call_gamma_data = {}
for T in expiry_periods:
    call_gamma_data[T] = []
    for K in K_range:
        gamma = calculate_call_gamma(S, K, r, T, sigma)
        call_gamma_data[T].append(gamma)

# Displaying Call Gamma values in a table
table_data = []
header_row = ["Strike Price / Time to Expiry"] + [f"{int(T*12)}M" for T in expiry_periods]
table_data.append(header_row)

for i, K in enumerate(K_range):
    row = [f"K={K:.2f}"]
    for T in expiry_periods:
        row.append(f"{call_gamma_data[T][i]:.6f}")
    table_data.append(row)

print(tabulate(table_data, headers="firstrow", tablefmt="pretty"))

# Plotting Call Gamma vs. Strike Price for different expiry periods
plt.figure(figsize=(10, 6))

for T in expiry_periods:
    plt.plot(K_range, call_gamma_data[T], marker='o', linestyle='-', label=f'{int(T*12)} Month Expiry')

plt.title('Call Gamma vs. Strike Price for Different Expiry Periods')
plt.xlabel('Strike Price')
plt.ylabel('Call Gamma')
plt.legend()
plt.grid(True)
plt.show()

Output

Gamma (Γ) for different Implied Volatility

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

# Function to calculate call gamma
def calculate_call_gamma(S, K, r, T, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    gamma = np.exp(-r * T) * norm.pdf(d1) / (S * sigma * np.sqrt(T))
    return gamma

# Parameters
S = 100  # Current stock price
K_range = np.arange(80, 121, 5)  # Strike price range at intervals of 5
r = 0.05  # Risk-free rate
expiry_period = 1/12  # Expiry period in years (1 month)
high_volatility = 0.5  # High volatility value
low_volatility = 0.3  # Low volatility value

# Calculate call gamma for strike prices at intervals of 5 and different volatility levels
call_gamma_high_vol = [calculate_call_gamma(S, K, r, expiry_period, high_volatility) for K in K_range]
call_gamma_low_vol = [calculate_call_gamma(S, K, r, expiry_period, low_volatility) for K in K_range]

# Plotting Call Gamma vs. Strike Price for high and low volatility
plt.figure(figsize=(8, 6))

plt.plot(K_range, call_gamma_high_vol, marker='o', linestyle='-', label=f'High Volatility ({high_volatility})')
plt.plot(K_range, call_gamma_low_vol, marker='o', linestyle='-', label=f'Low Volatility ({low_volatility})')

plt.title('Call Gamma vs. Strike Price for Different Volatility Levels')
plt.xlabel('Strike Price')
plt.ylabel('Call Gamma')
plt.legend()
plt.grid(True)
plt.show()

# Displaying Call Gamma values in a table for high and low volatility levels
table_data = []
for i, K in enumerate(K_range):
    table_data.append([f"K={K:.2f}", f"{call_gamma_high_vol[i]:.6f}", f"{call_gamma_low_vol[i]:.6f}"])

headers = ['Strike Price', f'High Volatility ({high_volatility})', f'Low Volatility ({low_volatility})']
print(tabulate(table_data, headers=headers, tablefmt="pretty"))

Theta (Θ) with Time to Expiry

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

# Function to calculate call theta
def calculate_call_theta(S, K, r, T, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    theta = (-S * norm.pdf(d1) * sigma / (2 * np.sqrt(T))) - r * K * np.exp(-r * T) * norm.cdf(d2)
    return theta

# Parameters
S = 100  # Current stock price
T_range = np.arange(0.25, 2.25, 0.25)  # Time to expiry range at intervals of 0.25
K = 100  # Strike price
r = 0.05  # Risk-free rate
volatility = 0.3  # Volatility value

# Calculate call theta for different time periods
theta_values = [calculate_call_theta(S, K, r, T, volatility) for T in T_range]

# Creating a list of tuples (Time to Expiry, Theta) for sorting
theta_tuples = list(zip(T_range, theta_values))
sorted_theta_tuples = sorted(theta_tuples, key=lambda x: x[1], reverse=True)  # Sorting based on Theta values in descending order

# Displaying Theta values in a table for different time periods
table_data = []
for T, theta in sorted_theta_tuples:
    table_data.append([f"Time={T:.2f}", f"{theta:.6f}"])

headers = ['Time to Expiry', 'Theta']
print(tabulate(table_data, headers=headers, tablefmt="pretty"))

# Plotting Theta vs. Time to Expiry
plt.figure(figsize=(8, 6))
plt.plot(T_range, theta_values, marker='o', linestyle='-', color='blue')
plt.title('Theta vs. Time to Expiry')
plt.xlabel('Time to Expiry')
plt.ylabel('Theta')
plt.grid(True)
plt.show()

Output

Theta (Θ) for different Implied Volatility

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

# Function to calculate call theta
def calculate_call_theta(S, K, r, T, sigma):
    d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    call_theta = (-S * norm.pdf(d1) * sigma / (2 * np.sqrt(T))) - r * K * np.exp(-r * T) * norm.cdf(d2)
    return call_theta

# Parameters
S = 100  # Current stock price
K_range = np.arange(80, 121, 5)  # Strike price range with intervals of 5
r = 0.05  # Risk-free rate
T = 1  # Time to expiry
high_volatility = 0.40  # High volatility value
low_volatility = 0.25  # Low volatility value

# Calculate call theta for high and low volatility with a strike price interval of 5
call_theta_high_vol = []
call_theta_low_vol = []

for K in K_range:
    call_theta_high = calculate_call_theta(S, K, r, T, high_volatility)
    call_theta_low = calculate_call_theta(S, K, r, T, low_volatility)
    call_theta_high_vol.append(call_theta_high)
    call_theta_low_vol.append(call_theta_low)

# Reverse the order of Theta values
call_theta_high_vol.reverse()
call_theta_low_vol.reverse()

# Plotting Call Theta vs. Strike Price for high and low volatility
plt.figure(figsize=(8, 6))
plt.plot(K_range, call_theta_high_vol, label='Call Theta (High Vol)', linestyle='-', color='blue')
plt.plot(K_range, call_theta_low_vol, label='Call Theta (Low Vol)', linestyle='--', color='red')

plt.title('Call Theta vs. Strike Price for Different Volatility Levels')
plt.xlabel('Strike Price')
plt.ylabel('Call Theta')
plt.legend()
plt.grid(True)
plt.gca().invert_yaxis()  # Invert the y-axis to display Theta in reverse order
plt.show()

# Displaying Call Theta values in a table for high and low volatility
table_data = []
for i, K in enumerate(K_range):
    table_data.append([f"K={K}", f"{call_theta_high_vol[i]:.6f}", f"{call_theta_low_vol[i]:.6f}"])

headers = ['Strike Price', 'Call Theta (High Vol)', 'Call Theta (Low Vol)']
print(tabulate(table_data, headers=headers, tablefmt="pretty"))