Babaoglu, Christoffersen, Heston y Jacobs (2018) introdujo un modelo GARCH de componentes con innovaciones gaussianas inversas y un núcleo de precios exponencialmente cuadrático allá por 2018. El artículo no debería estar detrás de un muro de pago, pero puedo intentar encontrarlo en otro lugar si lo está.
La fijación de precios en este modelo es muy conveniente. Al igual que en Heston y Nandi (2000), el modelo admite una función generadora de momentos exponencialmente afín, por lo que tenemos una fórmula casi analítica. Como referencia, las páginas relevantes para nosotros aquí son:
- Páginas 193-194: Dan la fórmula de precios, junto con la recursividad habitual (renuncio a escribirlas aquí porque son muy largas y no quiero arriesgarme a equivocarme);
- Página 202: Tiene un conjunto de parámetros para el modelo IG-GARCH(C) que puede utilizar para comprobar si la función tiene sentido.
Ahora, codifiqué la fórmula de precios en Python, pero no parece comportarse correctamente. En concreto, tengo funciones que sé que funcionan correctamente para el modelo de Heston y Nandi (2000), por ejemplo. Las obtuve del sitio web de Christoffersen, las traduje a Python y hago coincidir la salida de MATLAB hasta algo así como el noveno o décimo decimal. Así que elegí un ejemplo:
- tipo libre de riesgo = 0,05/365
- precio actual del subyacente = 100
- precio de ejercicio = 100
- volatilidad inicial = 0,21^2/365
- días hasta el vencimiento = 60 días
Para Heston y Nandi (2000), esa opción vale \$3.7778; For Black-Scholes-Merton, it's worth \$ 3.3968.
Así que espero algo en esas aguas. Pero, utilizando los parámetros para el proceso de riesgo neutral de la página 202, e inicializando ambos procesos de volatilidad ( $q_t$ y $h_t$ en el modelo) a 0,21^2/365, obtengo por encima de 47$, lo que es un sinsentido evidente. Podría ser un error de codificación, pero también podría ser algo estúpido sobre el escalado que no veo o podría ser que introdujera los parámetros incorrectos... Sigo revisando, pero tengo la nariz pegada y no veo dónde está el problema.
I eligió para publicar esto aquí porque necesito a alguien que esté lo suficientemente familiarizado con la fijación de precios de opciones para ver lo que está mal, si no es sólo un error de codificación. Por el lado bueno, si arreglamos el código de ejemplo aquí, todos en el foro podrán disfrutar de un código abierto para un modelo de última generación con características muy interesantes. Ten en cuenta que comenté otra forma de hacer el precio usando sólo una integral. De momento no importa. Mi código en Python:
import numpy as np
from numpy import sqrt, exp, log
from scipy.integrate import quad
# BCHJ2018, p.20
# mu_t,wq,rho1,ah,ch,rho2,aq,cq,eta
param = [-0.5, 2.415e-6, 0.745, 1.033e6, 9.682e-7,
0.989, 4.911e7, 4.660e-6, -5.399e-4]
# To try the functions
BSvol = 0.21
qt = BSvol**2/365
ht = BSvol**2/365
St = 100
K = 100
tau = 60
rF = 0.05/365
#==========================================================================#
def CF_IG_GARCH_C(u,St,rF,ht,qt,tau,param):
'''
Author: Stephane Surprenant, UQAM
Creation: 14/03/2020
Description: This function provides the generating function used in the
valuation of European call options for the IG-GARCH(C) model of Babaoglu,
Christoffersen, Heston and Jacobs (2018).
INPUTS DESCRIPTION
u : (float) Value over which the CGF is integrated
St: (float) Stock/index level at time t
ht: (float) Daily variance in the 2nd period (t+1)
(Vol.daily = Vol.yearly^2/365)
qt: (float) Daily long term variance in the 2nd period (t+1)
tau: (int) Time to maturity (days)
r : (float) Daily risk-free rate (rf.daily = rf.yearly/365)
param: (float) Array: [mu_t,wq,rho1,ah,ch,rho2,aq,cq,eta]
Note: Those are the risk-neutral component model risk-neutral parameters.
References: See Babaoglu, Christoffersen, Heston and Jacobs (2018).
REQUIRES: numpy (import sqrt, log, exp)
'''
# Assign parameter values
mu_t,wq,rho1,ah,ch,rho2,aq,cq,eta = param
mu = mu_t - eta**(-1) # (p.11, top)
# Complex argument
#u1 = u*1j
u1 = u
T = tau
# Matrices for the recursion (impose A(T)=B(T)=0)
Amat = np.zeros(shape=(T), dtype=complex)
Bmat = np.zeros(shape=(T), dtype=complex)
Cmat = np.zeros(shape=(T), dtype=complex)
e2 = eta**2
e4 = eta**4
# Initialize matrices at T-1
Amat[0] = u1*rF
Bmat[0] = mu*u1 + eta**(-2) - eta**(-2)*sqrt(1-2*eta*u1)
Cmat[0] = 0
# Recursion backward in time (first is last in the matrix)
for tt in range(1,T):
Amat[tt] = Amat[tt-1] + u1*rF + (wq - ah*e4 - aq*e4)*Bmat[tt-1] \
+ (wq - aq*e4)*Cmat[tt-1] \
-0.5*log(1 - 2*(ah+aq)*e4*Bmat[tt-1] \
- 2*aq*e4*Cmat[tt-1])
Bmat[tt] = u1*mu + (rho1 - (ch+cq)*eta**(-2) -(ah+aq)*e2)*Bmat[tt-1] \
- (cq*eta**(-2) + aq*eta**2)*Cmat[tt-1] + eta**(-2) \
- eta**(-2)*sqrt((1 - 2*(aq+ah)*e4*Bmat[tt-1]\
- 2*aq*e4*Cmat[tt-1])\
*(1 - 2*eta*u1-2*(cq+ch)*Bmat[tt-1] \
- 2*cq*Cmat[tt-1]))
Cmat[tt] = (rho2-rho1)*Bmat[tt-1] + rho2*Cmat[tt-1]
# g_t(u1,T) : (St**u1)*exp(A(t)+B(t)*h(t+1)+C(t)*q(t+1))
gt = exp(log(St)*u1 + Amat[tau-1] + Bmat[tau-1]*ht + Cmat[tau-1]*qt)
return(gt)
#==========================================================================#
def Price_IG_GARCH_C(St,K,rF,ht,qt,tau,param):
'''
Author: Stephane Surprenant, UQAM
Creation: 15/03/2020
Description: Valuation of European call options for the IG-GARCH(C) model
of Babaoglu, Christoffersen, Heston and Jacobs (2018) using IFT.
INPUTS DESCRIPTION
K : (float) Strike price
St: (float) Stock/index level at time t
ht: (float) Daily variance in the 2nd period (t+1)
(Vol.daily = Vol.yearly^2/365)
qt: (float) Daily long term variance in the 2nd period (t+1)
tau: (int) Time to maturity (days)
r : (float) Daily risk-free rate (rf.daily = rf.yearly/365)
param: (float) Array: [mu_t,wq,rho1,ah,ch,rho2,aq,cq,eta]
Note: Those are the risk-neutral component model risk-neutral parameters.
References: See Babaoglu, Christoffersen, Heston and Jacobs (2018).
REQUIRES: numpy (import sqrt, log, exp), scipy.integrate (quad)
'''
# Integrands
f1 = lambda u: np.real(K**(-u*1j)*\
CF_IG_GARCH_C(u*1j+1,St,rF,ht,qt,tau,param)/(St*u*1j))
f2 = lambda u: np.real(K**(-u*1j)*\
CF_IG_GARCH_C(u*1j,St,rF,ht,qt,tau,param)/(u*1j))
# Pricing formula (p.11)
cPrice = St*(0.5 + exp(-rF*tau)/np.pi*quad(f1,0,10000)[0]) \
- K*exp(-rF*tau)*(0.5+1/np.pi*quad(f2,0,10000)[0])
# =============================================================================
# t_Hk = lambda u: np.imag(CF_IG_GARCH(u-1j,St,ht,tau,rF,param) \
# *exp(-1j*u*log(K))/(1j*u+1))/u
# cPrice = 0.5*St + exp(-rF*tau)/np.pi*quad(t_Hk, 0, 1000)[0]
# =============================================================================
return(cPrice)
