Расчет и Визуализация Поверхности Подразумеваемой Волатильности (IV Surface) Опционов в Python Matplotlib

import numpy as np
from scipy.stats import norm
from scipy.optimize import brentq
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# --- 1. Модель Блэка-Шоулза-Мертона для колл-опциона ---
def black_scholes_call(S, K, T, r, sigma):
    """
    Рассчитывает цену европейского колл-опциона по модели Блэка-Шоулза-Мертона.

    Параметры:
    S (float): Текущая цена базового актива.
    K (float): Цена страйка опциона.
    T (float): Время до экспирации (в годах).
    r (float): Безрисковая процентная ставка (годовая).
    sigma (float): Подразумеваемая волатильность (годовая).

    Возвращает:
    float: Цена колл-опциона.
    """
    if T <= 0: # Опционы с нулевым или отрицательным временем до экспирации
        return max(0, S - K)
    
    d1 = (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    
    call_price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
    return call_price

# --- 2. Функция для расчета подразумеваемой волатильности (IV) ---
def implied_volatility_call(option_price, S, K, T, r, initial_vol=0.5, tol=1e-6, max_iter=100):
    """
    Рассчитывает подразумеваемую волатильность для европейского колл-опциона.

    Параметры:
    option_price (float): Рыночная цена опциона.
    S (float): Текущая цена базового актива.
    K (float): Цена страйка опциона.
    T (float): Время до экспирации (в годах).
    r (float): Безрисковая процентная ставка (годовая).
    initial_vol (float): Начальное предположение для волатильности.
    tol (float): Допуск для сходимости.
    max_iter (int): Максимальное количество итераций.

    Возвращает:
    float: Подразумеваемая волатильность или np.nan, если не найдена.
    """
    if T <= 0: # Для опционов с нулевым временем до экспирации IV не определена
        return np.nan

    # Функция, которую нужно обнулить для нахождения IV
    def price_diff(sigma):
        return black_scholes_call(S, K, T, r, sigma) - option_price

    # Границы для поиска волатильности
    low_vol = 1e-5  # Очень низкая волатильность
    high_vol = 5.0  # Очень высокая волатильность (500%)

    # Проверка на возможность существования IV
    # Если рыночная цена ниже минимальной (при sigma=0) или выше максимальной
    # (при очень высокой sigma), IV может быть вне диапазона или не существовать.
    # В таких случаях brentq может не найти корень или выдать ошибку.
    # Мы проверяем, что функция price_diff меняет знак в диапазоне [low_vol, high_vol].
    if price_diff(low_vol) * price_diff(high_vol) > 0:
        return np.nan

    try:
        # Используем метод Брента для нахождения корня
        iv = brentq(price_diff, low_vol, high_vol, xtol=tol, maxiter=max_iter)
        return iv
    except ValueError:
        return np.nan

# --- 3. Генерация синтетических данных опционов и расчет IV Surface ---
def generate_and_plot_iv_surface(S0=100, r=0.01, num_strikes=20, num_maturities=15):
    """
    Генерирует синтетические данные опционов, рассчитывает IV и строит 3D-поверхность.

    Параметры:
    S0 (float): Текущая цена базового актива.
    r (float): Безрисковая процентная ставка.
    num_strikes (int): Количество страйков для генерации.
    num_maturities (int): Количество сроков до экспирации для генерации.
    """
    print(f"Генерация данных для S0={S0}, r={r}")

    # Диапазоны страйков и сроков до экспирации
    strikes = np.linspace(S0 * 0.8, S0 * 1.2, num_strikes) # От 80% до 120% от S0
    maturities = np.linspace(0.1, 2.0, num_maturities) # От 0.1 года до 2 лет

    # Создаем сетку для страйков и сроков
    Strikes_mesh, Maturities_mesh = np.meshgrid(strikes, maturities)

    # Инициализируем массив для хранения подразумеваемой волатильности
    iv_surface = np.zeros_like(Strikes_mesh, dtype=float)

    # --- Симуляция рыночных цен опционов с реалистичным IV "смещением" ---
    # Создадим базовую IV поверхность с улыбкой/перекосом
    base_iv = 0.20 # Базовая волатильность
    skew_factor = 0.005 # Фактор перекоса (skew)
    smile_factor = 0.0005 # Фактор улыбки (smile)

    # Добавляем шум для реалистичности
    np.random.seed(42) # Для воспроизводимости
    noise = np.random.normal(0, 0.01, iv_surface.shape)

    print("Расчет подразумеваемой волатильности...")
    for i in range(num_maturities):
        for j in range(num_strikes):
            K = Strikes_mesh[i, j]
            T = Maturities_mesh[i, j]

            # Моделируем "истинную" IV с учетом страйка и срока
            # Чем дальше от ATM, тем выше IV (улыбка)
            # Чем ниже страйк, тем выше IV (перекос)
            simulated_iv = base_iv + \
                           skew_factor * (S0 - K) + \
                           smile_factor * (K - S0)**2 / S0**2 + \
                           0.02 * (1 - np.exp(-T/0.5)) # Зависимость от срока
            
            simulated_iv = max(0.05, simulated_iv) # Минимальная IV
            simulated_iv = min(1.0, simulated_iv) # Максимальная IV

            # Рассчитываем "рыночную" цену опциона на основе этой симулированной IV
            # Добавляем небольшой шум к цене, чтобы имитировать реальные данные
            simulated_option_price = black_scholes_call(S0, K, T, r, simulated_iv) * (1 + noise[i, j])

            # Рассчитываем IV из симулированной рыночной цены
            iv = implied_volatility_call(simulated_option_price, S0, K, T, r)
            iv_surface[i, j] = iv
    
    # Заменяем NaN значения на интерполированные или средние, чтобы не было "дыр" на графике
    # В реальных данных это может потребовать более сложной обработки
    if np.any(np.isnan(iv_surface)):
        print("Внимание: Обнаружены NaN значения в IV поверхности. Попробуйте скорректировать параметры или диапазон.")
        # Простая замена NaN на среднее значение (для демонстрации)
        mean_iv = np.nanmean(iv_surface)
        iv_surface[np.isnan(iv_surface)] = mean_iv

    # --- 4. Визуализация IV Surface ---
    fig = plt.figure(figsize=(12, 9))
    ax = fig.add_subplot(111, projection='3d')

    # Построение поверхности
    surf = ax.plot_surface(Strikes_mesh, Maturities_mesh, iv_surface, cmap='viridis')

    # Настройка осей и заголовка
    ax.set_xlabel('Страйк (K)')
    ax.set_ylabel('Время до экспирации (T, годы)')
    ax.set_zlabel('Подразумеваемая волатильность (IV)')
    ax.set_title(f'Поверхность подразумеваемой волатильности (IV Surface) для S0={S0}')

    # Добавление цветовой шкалы
    fig.colorbar(surf, shrink=0.5, aspect=5, label='IV')

    # Настройка угла обзора для лучшей видимости
    ax.view_init(elev=20, azim=-120) # elev - угол возвышения, azim - азимут

    plt.tight_layout()
    plt.show()

# --- Запуск генерации и построения --- 
if __name__ == '__main__':
    generate_and_plot_iv_surface()