Problema con la fijación del precio de una opción de compra utilizando el modelo de Monte Carlo Vasicek

Question

Estoy intentando fijar el precio de una opción de compra sobre un cupón cero según el modelo Vasicek utilizando el método de Monte Carlo:

$$C_0 = B(0,\theta) \ \mathbb{E}^{\mathbb{Q}_T}[(B(\theta,T)-K)^{+}]$$

El problema es que el código que he escrito para esto está dando un precio demasiado diferente de la respuesta esperada que es 0,095 basado en las fórmulas cerradas.

¿Qué puedo estar haciendo mal?

Los supuestos subyacentes y las ecuaciones que estoy utilizando:

Sea $(, (\mathcal{F}_t)_t, \mathbb{Q})$ tal que $\mathbb{Q}$ es la probabilidad neutral al riesgo.

Modelo Vasicek: $dr_t = a(b - r_t ) dt + \sigma dW_t^{\mathbb{Q}}$

Precio del bono cupón cero: $B(T,S) = \tau(T,S) \ exp(-\eta(T,S)\ r_t)$

Desde $r_t$ es gaussiano bajo $\mathbb{Q}_t$ (medida de avance): $$r(t) \sim \mathcal{N}(r_0e^{-at} + b(1-e^{-at}) + \frac{\sigma(1 - e^{-at})}{a^2} - \frac{\sigma(1 - e^{-2at})}{2a^2} \ , \ \sigma^{2} \frac{1 - e^{-2at}}{2a}) $$

Y mi código basado en lo anterior: (uno puede simplemente copiar y pegar y luego ejecutar el código de abajo asumiendo que las librerías externas de Python requeridas ya están instaladas).

import numpy as np
import pylab as plt
from random import gauss
from math import exp, sqrt, log
from scipy.stats import norm
a = 0.2
b = 0.1
sigma = 0.02
r = 0.1
t = 0
theta = 0.5
T = 1
K = 0.9

#Distribution of rt
r_0 = r
m_rt = r * exp(-a*theta) + b * (1 - exp(-a*theta)) + sigma * ((1 - exp(-a*theta)) / a**2 ) - sigma * ((1 - exp(-2*a*theta))/2*(a**2))
v_rt = (sigma**2) * (1 - exp(2*a*theta))/ 2*a

eta = (1 - exp(-a*(T-theta)))/a
tau = np.exp( (b - (sigma**2)/(2*(a**2))) * (eta - T + theta) - (0.25 * (sigma**2) * (eta**2) / a) )

def VasicekZCPrice(t,T,a,b,sigma,r):
    B = ( 1- exp( -a*(T-t) ) )/ a
    A = exp( ( (b-(0.5*sigma**2)) / a**2 )*( B - T + t ) - ( 0.25*sigma**2 / a ) * B**2)
    ZC = A*exp(-B*r)
    return ZC

def monte_carlo_call(theta, T, a, b, sigma, r, K, nbSimul = 10000, steps = 100):
    dt = (T-theta)/steps
    r__T = r * np.exp(-a*dt) +  b* (1-np.exp(-a*dt)) + (sigma*(1-exp(-a*dt))/a**2) - sigma * ((1-np.exp(-2*a*dt))/2*a**2) + sigma*np.sqrt((1-np.exp(-2*a*dt))/2*a) * np.random.normal(size=(steps,nbSimul))
    B__T = tau * np.exp(-eta*r__T)
    payoffs = np.maximum( B__T - K, 0 )
    option_price = VasicekZCPrice(t, T,a,b,sigma,r) * np.mean(payoffs)
    return option_price

monte_carlo_call(theta, T, a, b, sigma, r,K ,  nbSimul = 10000, steps = 100)

que me da:

0.03686737433679907

Se puede acceder a una versión ampliada del código y a más explicaciones a través de ce lien en Google Colab.

¡Muchas gracias!

Answer 1

Para asegurarme de que entiendo el problema: usted está tratando de fijar el precio de una opción de compra que vence en el momento 0,5, que se ejercerá en un bono de cupón cero nocional unitario con un vencimiento de 1,0 a un strike (precio) de 0,9.

Seguiré la notación de "Interest Rate Models - Theory and Practice" de Brigo & Mercurio (2006) que es ligeramente diferente a su notación.

De la página 58 del libro mencionado, obtenemos la dinámica del modelo Vasicek como:

$$ dr(t)=k[\theta-r(t)]dt+\sigma dW(t) $$

Y nosotros elegimos los parámetros:

k = 0.2
theta = 0.1
sigma = 0.02
r_0 = 0.1

time_now = 0
expiry = 0.5
maturity = 1
strike = 0.9

En QuantLib (C++), la función de interés ya está implementada. Por desgracia, esta función no está expuesta en la vinculación con Python, pero podemos utilizarla como referencia:

Vasicek vasicek(0.1,0.2,0.1,0.02);
std::cout << std::setprecision(16);
std::cout << vasicek.discountBondOption(Option::Type::Call,0.9,0.5,1.0) << std::endl;

Devuelve un valor de 0.0487763759011719 . Veamos si podemos obtener este valor con Python.

A partir de la página 59 del mismo libro anterior, podemos calcular el precio de un bono cupón cero utilizando:

$$ P(t,T)=A(t,T)e^{-B(t,T)r(t)} $$ donde $$ A(t,T)=\exp\left\{\left(\theta-\frac{\sigma^{2}}{2k^{2}}\right)\left[B(t,T)-T+t\right]-\frac{\sigma^{2}}{4k}B(t,T)^{2}\right\} $$ et $$ B(t,T)=\frac{1}{k}\left[1-e^{-k(T-t)}\right] $$

Para calcular estas funciones en Python, importamos los siguientes paquetes:

import numpy as np
from math import exp, sqrt, log
from scipy.stats import norm

e implementar las funciones de la siguiente manera:

def B(t,T):
    return (1/k)*(1-exp(-k*(T-t)))
def A(t,T):
    return exp((theta-pow(sigma,2)/(2*pow(k,2)))*(B(t,T)-T+t)-(pow(sigma,2)/(4*k))*pow(B(t,T),2))
def P(t,T,r_t):
    return A(t,T)*exp(-B(t,T)*r_t)

De la página 60, tenemos la fórmula de la opción de compra de bonos cupón cero:

$$ ZBO(t,T,S,X)=P(t,S)\Phi(h)-XP(t,T)\Phi(h-\sigma_{p}) $$

donde

$$ \sigma_{p}=\sigma\sqrt{\frac{1-e^{-2k(T-t)}}{2k}}B(T,S) $$

$$ h=\frac{1}{\sigma_{p}}\ln\frac{P(t,S)}{P(t,T)X}+\frac{\sigma_{p}}{2} $$

que implementamos como:

def ZBO(t,T,S,X,r_t):
    sigma_p=sigma*sqrt((1-exp(-2*k*(T-t)))/(2*k))*B(T,S)
    h=(1/sigma_p)*log(P(t,S,r_t)/(P(t,T,r_t)*X))+sigma_p/2
    return P(t,S,r_t)*norm.cdf(h, 0, 1)-X*P(t,T,r_t)*norm.cdf(h-sigma_p, 0, 1)

El valor de la opción calculado mediante ZBO(time_now,expiry,maturity,strike,r_0) sale como 0.0487763759011719 .

Este valor es el mismo que el calculado con QuantLib. No estoy seguro de dónde sacaste tu fórmula de forma cerrada, pero parece un poco apagado.

La realización de una simulación Monte Carlo es un poco complicada porque no sólo nos interesa el resultado simulado, sino también el resultado final. $r(t)$ de $0$ a $T$ sino también el tipo de interés integrado $\int_0^Tr(t)dt$ . Esto se debe a que estamos tratando de evaluar:

$$ e^{-\int_0^Tr(t)dt}\left[P(T,S)-X\right]^+ $$

Mientras que el valor del bono cupón cero entre paréntesis sólo depende del tipo de interés terminal $r(T)$ El factor de descuento inicial depende del tipo de interés al contado instantáneo, es decir, del tipo de interés a corto plazo.

Consideremos ahora la siguiente aproximación:

$$ \int_0^Tr(t)dt\approx\sum_{i=0}^{n}r(i\cdot \Delta t)\cdot \Delta t $$

Esto significa que dividiremos el tiempo en $n$ piezas y utilizar el tipo de interés al principio de cada período discreto para aproximar el tipo integrado.

Lo primero y más importante es simular la propia tasa. A partir de la página 59, tenemos la distribución del tipo a corto bajo la medida de riesgo neutro ( $\mathbb{Q}$ ) :

$$ \mathbb{E}\left[ r(t) \middle| \mathcal{F}_{s} \right] =r(s)e^{-k(t-s)}+\theta\left(1-e^{-k(t-s)}\right) $$

$$ \text{Var}\left[ r(t) \middle| \mathcal{F}_{s} \right] =\frac{\sigma^{2}}{2k}\left[1-e^{-2k(t-s)}\right] $$

cuya implementación en Python es la siguiente:

def E_r(s,t,r_s):
    return r_s*exp(-k*(t-s))+theta*(1-exp(-k*(t-s)))
def Var_r(s,t,r_s):
    return (pow(sigma,2)/(2*k))*(1-exp(-2*k*(t-s)))

En la simulación Monte Carlo, utilizamos pasos equidistantes de tamaño $\Delta t$ . Así, uno de los argumentos que utilizamos para llamar a estas funciones es 0,dt . El tipo de interés introducido en la fórmula se actualiza en cada paso. Esto nos lleva a la siguiente implementación de la simulación Monte Carlo:

np.random.seed(2)
n_steps=100
n_sim=10000

def mc_single(t,T,S,X):
    r=r_0
    dt=(T-t)/n_steps
    time=t
    interest_rate_integral=0
    for x in range(n_steps):
        interest_rate_integral=interest_rate_integral+r*dt
        r=np.random.normal(E_r(0.0,dt,r),sqrt(Var_r(0.0,dt,r)), 1)
        time=time+dt
    disc=exp(-interest_rate_integral)
    return disc*max(P(T,S,r)-X,0)

def mc(t,T,S,X):
    sum=0
    for x in range(n_sim):
        sum=sum+mc_single(t,T,S,X)
    return sum/n_sim

Estoy poniendo la semilla para que sea posible reproducir el resultado. El valor simulado calculado mediante mc(time_now,expiry,maturity,strike) sale como 0.0488357730948401 .

No es exactamente igual, pero se aproxima bastante teniendo en cuenta la simplicidad de este algoritmo.