Babaoglu, Christoffersen, Heston and Jacobs (2018) introduced a component GARCH model with inverse Gaussian innovations and an exponentially quadratic pricing kernel back in 2018. The article shouldn't be behind a paywall, but I can try to find it elsewhere if it is.
Pricing in this model is very convenient. As in Heston and Nandi (2000), the model admits an exponentially affine moment generating function, so we have a quasi-analytical formula. For reference, the relevant pages for us here are:
- Pages 193-194: They give the pricing formula, along with the usual recursion (I forgo writing them here because they are very long and I don't want to risk making a mistake);
- Page 202: You have a set of parameters for the IG-GARCH(C) model that you can use to check if the function makes sense.
Now, I coded the pricing formula in Python, but it doesn't seem to behave properly. Specifically, I have functions that I know work properly for the Heston and Nandi (2000) model, for example. I got them from Christoffersen's website, translated them in to Python and I match the MATLAB output up to something like the 9th or 10th decimal. So, I picked an example:
- risk free rate = 0.05/365
- current price of the underlying = 100
- strike price = 100
- initial volatility = 0.21^2/365
- days to maturity = 60 days
For Heston and Nandi (2000), that option is worth \$3.7778; For Black-Scholes-Merton, it's worth \$3.3968.
So, I'm expecting something in those waters. But, using the parameters for the risk-neutral process from page 202, and initializing both volatility processes ($q_t$ and $h_t$ in the model) at 0.21^2/365, I get above 47$ which is obvious nonsense. It could be a coding mistake, but it could also be something stupid about scaling that I don't see or it could be that I inputed the wrong parameters... I keep checking, but my nose is stuck in it and I just don't see where lies the problem.
I chose to post this here because I need someone who is familiar enough with option pricing to see what is wrong, if it isn't just a coding mistake. On the bright side, if we fix the sample code here, everyone on the forum will get to enjoy open source code for a state-of-the-art model with very cool features. Note that I commented another way to do the pricing using just one integral. It doesn't really matter for the moment. My Python code:
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)
