Comprensión de Monte Carlo para resolver el precio de las opciones con volatilidad local

Question

He leído esta pregunta fijación de precios mediante el modelo de volatilidad local dupire que parece tener una respuesta desde aquí https://www.csie.ntu.edu.tw/~d00922011/python/cases/LocalVol/DUPIRE_FORMULA.PDF

Ambos dicen simplemente que hay que utilizar métodos de monte carlo o de diferencias finitas una vez que tenemos la volatilidad local para valorar las opciones, después de discutir un poco más a fondo cómo encontrar la volatilidad local en sí.

Estoy buscando una descripción más granular de este último paso que es encontrar el precio de la opción en sí.

Mi intento

En el punto en el que estoy atascado, tengo todos los valores necesarios para resolver la ecuación dúplex: $dS_t=\mu(S_t,t) S_tdt+\sigma(S_t;t) S dW$ donde $S_t|_{t=0} = S_0$ . Supongamos que tengo bien definidos los términos de deriva y volatilidad local y que conozco también el precio al contado. A partir de aquí, puedo utilizar una solución monte carlo para encontrar $S_t$ a todos los conocidos $t$ . Sin embargo, $S_t$ no es el precio de la opción, por lo que me pregunto cómo ha podido encontrar estos precios la persona que ha escrito el artículo/ha formulado la pregunta vinculada.

Desde mi (limitada) comprensión de los métodos monte carlo y las ecuaciones diferenciales estocásticas, $S_t$ encontrada mediante los métodos MC es una realización de la variable aleatoria que es el precio subyacente de esta ecuación.

Sé que si encuentro la distribución de esta variable aleatoria, $p(x)$ En cada paso de tiempo podría ponerle precio a mi opción en cada paso de tiempo $\big($$ C(K,t) = E[(S_{t}-K)^+] = \int_{K}^infty p(x)(S_{t}-K)dx $$\big)$ .

He leído algunas descripciones sobre cómo resolver las integrales de Ito y Stratonovich utilizando métodos de monte carlo, pero éstas simplemente me dicen cómo evaluar estas integrales, lo cual ya sé hacer; quiero encontrar la probabilidad. ¿Cuál es la mejor manera de hacerlo? ¿Existen mejores métodos o paquetes de python que deba utilizar en su lugar?

Answer 1

Voy a mostrar un ejemplo codificado para mostrar cómo se puede intentar esto, utilizando el puerto de python del QuantLib biblioteca. Todo parecerá un poco mecánico, pero espero que sea instructivo. Hay un poco de código de configuración requerida (especificando los tipos de interés, spot, etc., y también una función de utilidad que utilizo para trazar las superficies), he empujado todo esto a la parte inferior de mi respuesta, pero lo necesitará en la parte superior de su script/notebook.

Para construir un modelo de vol local, voy a necesitar una superficie de vol implícita de los datos del mercado. No sé qué datos estás utilizando, así que como ejemplo muy trivial, voy a simular una secuencia de sonrisas SABR a diferentes vencimientos, y asumir que es la superficie de vol implícita del mercado. El código para hacer eso es:

strikes = [70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0]
expirations = [ql.Date(1, 7, 2021), ql.Date(1, 9, 2021), ql.Date(1, 12, 2021), ql.Date(1, 6, 2022)]
vol_matrix = ql.Matrix(len(strikes), len(expirations))

# params are sigma_0, beta, vol_vol, rho
sabr_params = [[0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6]]

for j, expiration in enumerate(expirations):
    for i, strike in enumerate(strikes):
        tte = day_count.yearFraction(today, expiration)
        vol_matrix[i][j] = ql.sabrVolatility(strike, spot, tte, *sabr_params[j])

implied_surface = ql.BlackVarianceSurface(today, calendar, expirations, strikes, vol_matrix, day_count)
implied_surface.setInterpolation('bicubic')
vol_ts = ql.BlackVolTermStructureHandle(implied_surface)

process = ql.BlackScholesMertonProcess(ql.QuoteHandle(ql.SimpleQuote(spot)), dividend_ts, flat_ts, vol_ts)

plot_vol_surface(implied_surface, plot_years=np.arange(0.1, 1, 0.1))
plt.title("Implied Vol")

Esto también traza la superficie del vol para que podamos verlo:

Ahora voy a intentar fijar el precio de una opción vainilla a 6 meses de vencimiento, utilizaré un motor de fijación de precios por diferencias finitas y un motor de fijación de precios monte carlo. Primero, vamos a consultar la superficie y ver qué vol esperamos:

implied_surface.blackVol(0.5, 90)

devuelve 0,0787116605540102, es decir, un 7,87% de IV

En QuantLib, configuramos los motores de opciones y de precios de forma independiente. Primero, la opción:

strike = 90.0
option_type = ql.Option.Call

maturity = today + ql.Period(6, ql.Months)
europeanExercise = ql.EuropeanExercise(maturity)

payoff = ql.PlainVanillaPayoff(option_type, strike)

european_option = ql.VanillaOption(payoff, europeanExercise)

y ahora configuraré y fijaré el precio utilizando cada uno de los dos motores de fijación de precios (se puede encontrar un recurso inestimable sobre ellos en el ReadTheDocs página mantenida por uno de los otros carteles aquí).

En primer lugar, el FD pricer (tenga en cuenta que TENEMOS que activar el parámetro final aquí, que especifica el uso de Dupire local vol...):

tGrid, xGrid = 3000, 400
fd_engine = ql.FdBlackScholesVanillaEngine(process, tGrid, xGrid, 0, ql.FdmSchemeDesc.Douglas(), True)

european_option.setPricingEngine(fd_engine)
fd_price = european_option.NPV()
fd_price

el precio es de 7,699526511916783

Y ahora el MC pricer (para el MC en vol local, un número suficiente de pasos es vital, experimente con la reducción de esto y vea lo que sucede con el precio):

steps = 36
rng = "lowdiscrepancy" # could use "pseudorandom"
numPaths = 2**15

mc_engine = ql.MCEuropeanEngine(process, rng, steps, requiredSamples=numPaths)

european_option.setPricingEngine(mc_engine)
mc_price = european_option.NPV()
mc_price

el precio es de 7,684257656799339

Y como comparación final, calcularemos los vols que implica cada uno de estos precios, para comparar con el vol de la superficie anterior:

# Compare calculated volatility to vol from the surface
european_option.impliedVolatility(fd_price, process), european_option.impliedVolatility(mc_price, process)

(0,07872046080861067, 0,077079492498713), ambos bastante cercanos al vol que esperábamos.

Espero que esto sea útil. Hay cientos de pequeñas suposiciones que he hecho en la respuesta, sin embargo, que podrían no ser apropiadas para tu situación... hazme saber si quieres que cambie algo.

Código de caldera necesario para ejecutar los fragmentos anteriores:

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.mplot3d import Axes3D
import QuantLib as ql

calc_date = ql.Date(21, ql.December, 2020)

def plot_vol_surface(vol_surface, plot_years=np.arange(0.1, 3, 0.1), plot_strikes=np.arange(70, 130, 1), funct='blackVol'):
    if type(vol_surface) != list:
        surfaces = [vol_surface]
        functs = [funct]
    else:
        surfaces = vol_surface

        if type(funct) != list:
            functs = [funct] * len(surfaces)
        else:
            functs = funct

    fig = plt.figure(figsize=(10, 6))
    ax = fig.gca(projection='3d')
    X, Y = np.meshgrid(plot_strikes, plot_years)
    Z_array, Z_min, Z_max = [], 100, 0

    for surface, funct in zip(surfaces, functs):
        method_to_call = getattr(surface, funct)

        Z = np.array([method_to_call(float(y), float(x)) 
                      for xr, yr in zip(X, Y) 
                          for x, y in zip(xr, yr)]
                     ).reshape(len(X), len(X[0]))

        Z_array.append(Z)
        Z_min, Z_max = min(Z_min, Z.min()), max(Z_max, Z.max())

    # In case of multiple surfaces, need to find universal max and min first for colourmap
    for Z in Z_array:
        N = (Z - Z_min) / (Z_max - Z_min)  # normalize 0 -> 1 for the colormap
        surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0.1, facecolors=cm.coolwarm(N))

    m = cm.ScalarMappable(cmap=cm.coolwarm)
    m.set_array(Z)
    plt.colorbar(m, shrink=0.8, aspect=20)
    ax.view_init(30, 300)

# Simple World State
spot = 100
rate_dom = 0.02
rate_for = 0.05

today = ql.Date(1, 6, 2021)

calendar = ql.NullCalendar()
day_count = ql.Actual365Fixed()

# Set up some risk-free curves
riskFreeCurveDom = ql.FlatForward(today, rate_dom, day_count)
riskFreeCurveFor = ql.FlatForward(today, rate_for, day_count)

flat_ts = ql.YieldTermStructureHandle(riskFreeCurveDom)
dividend_ts = ql.YieldTermStructureHandle(riskFreeCurveFor)

Actualización - Interpolación de Andraesen-Huge (nb. He cambiado ligeramente la función gráfica de arriba):

strikes = np.linspace(70, 130, 10)
expirations = [ql.Date(1, 7, 2021), ql.Date(1, 8, 2021), ql.Date(1, 9, 2021), ql.Date(1, 12, 2021), ql.Date(1, 6, 2022)]

# params are sigma_0, beta, vol_vol, rho
sabr_params = [[0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6]]

# Now try Andraeson Huge calibration
calibration_set = ql.CalibrationSet()

for i, strike in enumerate(strikes):
    for j, expiration in enumerate(expirations):
        tte = day_count.yearFraction(today, expiration)
        payoff = ql.PlainVanillaPayoff(ql.Option.Call, strike)
        exercise = ql.EuropeanExercise(expiration)
        vol = ql.sabrVolatility(strike, spot, tte, *sabr_params[j])

        calibration_set.push_back((ql.VanillaOption(payoff, exercise), ql.SimpleQuote(vol)))

ah_interpolation = ql.AndreasenHugeVolatilityInterpl(calibration_set, \
                              ql.QuoteHandle(ql.SimpleQuote(spot)), flat_ts, dividend_ts)
ah_surface = ql.AndreasenHugeVolatilityAdapter(ah_interpolation)

plot_vol_surface([implied_surface, ah_surface], plot_years=np.arange(0.2, 1, 0.05), plot_strikes=np.arange(80, 120, 2))