Numerical Integration Grids for Heston–Nandi Option Pricing

# Pricing options by IFT under the Heston and Nandi (2000) model: odd behavior


# Pricing options by IFT under the Heston and Nandi (2000) model: odd behavior












I am working on option pricing using GARCH models and, currently, I am coding the pricing of options under the Heston and Nandi (2000) model. This model admits a quasi analytical formula for pricing options. As a benchmark, I originally worked out some code that mimicks the behavior of the functions provided by Peter Christoffersen's MATLAB toolbox. You can check it up on his website. I do not fall exactly on the same values he finds, but it's pretty darn close.

Now, for those of you who ever did this, you are aware that the coefficients in the characteristic function are computed through recursions going backward in time. Fortunately, for a given day and a given maturity, those coefficients are the same (because the risk-free rate is the same for all those contracts). So, I did a slightly modified version of my code where I pre-compute the coefficients. I simply define a grid of points for the integration and get the coefficients for all the points on the grid. Then, I can compute the characteristic function and use it for pricing through IFT.

The difference? Well, the original code uses a python lambda function with scipy.integrate.quad while the newer version uses a predefined grid phi with scipy.integrate.trapz. I also tried with scipy.integrate.simps and results are similar to those with trapz. In essence, my old functions match very well the MATLAB Toolbox, but my new functions are quite a bit off. Here is the code you need to reproduce the problem. Hopefully, someone will see the problem... I really don't know what is going on here.

```
import numpy as np
from   numpy import sqrt, exp, log
from   scipy.integrate import quad, trapz, simps
#%%

PHI = []
FIRST_PASS = True
LIMIT = 1000

def CF_HN(u,St,ht,r,tau,param):
    Description: This computes the characteristic function of the Heston and
    Nandi (2000) model.

    INPUTS       DESCRIPTION
    u :          (float) Value over which the CF 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)
    tau:         (int)   Time to maturity (days)
    r  :         (float) Daily risk-free rate (rf.daily = rf.yearly/365)

    param:       (array) [alpha, beta, gamma, h0, Lambda] where h0 is the
                          target variance.

    REQUIRES: numpy (import sqrt, log, exp)
    '''
    if FIRST_PASS:
       PHI.append(u)

    # Assign parameter values
    alpha, beta, gamma, h0, Lambda = param
    # Variance targeting imposes
    omega = (1 - beta - alpha*gamma**2)*h0 - alpha

    # Initialize
    u1 = u*1j
    T  = tau
    Amat = np.zeros(shape=(T), dtype=complex)
    Bmat = np.zeros(shape=(T), dtype=complex)

    # At time T-1
    Amat[0] = u1*r
    Bmat[0] = Lambda*u1 + 0.5*u1**2

    # Recursion backward in time (first is last in the matrix)
    for tt in range(1,T):
        Amat[tt] = Amat[tt-1] + u1*r + Bmat[tt-1]*omega \
                          - 0.5*log(1 - 2*alpha*Bmat[tt-1])
        Bmat[tt] = u1*(Lambda + gamma) - 0.5*gamma**2 + beta*Bmat[tt-1] \
               + (0.5*(u1-gamma)**2)/(1-2*alpha*Bmat[tt-1])

    Psi = exp(log(St)*u1 + Amat[T-1] + Bmat[T-1]*ht)
    return(Psi)

def Price_HN(St,K,ht,r,tau,param):
    '''
    Description: This function computes the price of European call options for
    the Heston and Nandi model 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)
    tau:         (int)   Time to maturity (days)
    r  :         (float) Daily risk-free rate (rf.daily = rf.yearly/365)

    param:       (array) [alpha, beta, gamma, h0, Lambda] where h0 is the
                         target variance.

    REQUIRES: numpy (import sqrt, log, exp), scipy.integrate (quad)
    '''

    t_Hk = lambda u: np.imag(CF_HN(u-1j,St,ht,r,tau,param)*exp(-1j*u*log(K))/(1j*u+1))/u
    cPrice = 0.5*St + exp(-r*tau)/np.pi*quad(t_Hk, 0, LIMIT)[0]

    return(cPrice)

#%% Modified versions of functions
if __name__ == '__main__':
# MATLAB CHECK:
St = np.arange(80,120.5,0.5)
K  = 100
r = 0.05/365
ht = 0.21**2/365
tau = 60

param = [7.83e-07, 0.8810, 378, 1.1100e-04, -0.5]
alpha, beta, gamma, lambda_ = np.hstack( (param[0:3], param[4]) )

sigma2      = param[3]
persistence = beta + alpha*gamma**2

param_ = np.array( [sigma2, alpha, persistence, gamma, lambda_] )

htp1 = 0.21**2/365
phi  = np.array( PHI )

def CharCoef(r, tau, param_, phi):

    sigma2, alpha, persistence, gamma, lambda_ = param_
    omega = (1 - beta - alpha*gamma**2)*sigma2 - alpha

    ui = (phi - 1j)*1j

    A = np.zeros(shape=(tau,len(ui)), dtype=complex)
    B = np.zeros(shape=(tau,len(ui)), dtype=complex)

    # At time DTM-1
    A[0] = ui*r
    B[0] = lambda_*ui + 0.5*ui**2
    # A[0] = 0*ui
    # B[0] = 0*ui

    # Recursion backward in time
    for dtm in range(1,tau):
        A[dtm] = A[dtm-1] + ui*r + B[dtm-1]*omega \
                          - 0.5*log( 1 - 2*alpha*B[dtm-1] )
        B[dtm] = ui*(lambda_ + gamma) - 0.5*gamma**2 +beta*B[dtm-1] \
                 + ( 0.5*(ui-gamma)**2 )/( 1-2*alpha*B[dtm-1] )

    return [A,B]

def CharFunc(St, tau, htp1, phi, A, B):
    ui = (phi - 1j)*1j

    Psi = exp(log(St)*ui + A[tau-1] + B[tau-1]*htp1)

    return Psi

def CPrice_HN(St, K, r, tau, Psi, phi):

    # integrand = np.imag( Psi*exp( -1j*phi*log(K) )/( 1j*phi+1 ) )/phi
    integrand = np.imag( Psi*exp(-1j*phi*log(K))/(1j*phi+1) )/phi
    CPrice    = 0.5*St + exp( -r*tau )/np.pi * trapz(integrand, phi)

    return CPrice
# np.imag(CF_HN(u-1j,St,ht,r,tau,param)*exp(-1j*u*log(K))/(1j*u+1))/u
A,B = CharCoef(r, tau, param_, phi)

p_old = np.array( [Price_HN(s,K,ht,r,tau,param) for s in St] )
p_new      = np.zeros( shape=p_old.shape )

count = 0
for s in St:
    Psi          = CharFunc(s, tau, htp1, phi, A, B)

    p_new[count] = CPrice_HN(s, K, r, tau, Psi, phi)
    count += 1

import pandas as pd

data = pd.DataFrame( {'Old Method': p_old, 'New Method': p_new } )
data.head()
```

Two notes: (1) Yes, it is bizarre to look at a range of prices, but it's intended to replicated what Christoffersen does in his Toolbox; (2) In the new version, I move back and forth between two sets of parameters. It's only for compatibility with other stuff I am doing.

EDIT: As it turns out, if I used a different integration grid, say

```
a    = np.arange(1e-20, 1, 0.001)
b    = np.arange(1,250,0.1)
phi_ = np.hstack( (a,b) )
```

I do recover nearly identical prices. I got this idea by noticing that the integrand has a much large absolute value near 0 and after about 100 or 200, it's worth almost nothing. So, I refined the grid and this gives me prices whose difference with old prices are on the order of 1e-7... In other words, it's negligible. This creates a new problem: how the hell do I justify sanely a choice for an integration grid that might look like this?