Correcting Risk-Neutral Drift in Inverse Gaussian Component GARCH Option Pricing

# Pricing options using the IG component GARCH model of BCHJ(2018)


# Pricing options using the IG component GARCH model of BCHJ(2018)












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)
```

## Answer by Stéphane (score 1, accepted)

https://quant.stackexchange.com/a/51808

I figured out the source of the problem today and it was indeed the stupid mistake that comes along with working on something like this late at night. In the paper, the repport risk-neutral estime is $\tilde{\mu} = \mu + 1/\eta$. Since $\eta$ is so damn small and negative, I was grossly understating the value of $\mu$...

I corrected the mistake directly above so that anyone interested in using the code may do it with ease. Now, as far as I can tell, the function seems to work adequately -- at least, if it doesn't, it spits out weirdly sensible prices when I compare them to other functions I know to work correctly.