0
$\begingroup$

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?

Improve this question
$\endgroup$
5
  • $\begingroup$ Essentially, the recursion is the one from Duffie/Pan/Singleton's 2000 Econometrica paper, no? I build a Matlab Toolbox for that, too. I remember that you need to take very good care around the Riccati integrals. I think I even added some Runge-Kutta methods to it just to stabilise these difference equations. Also, in your first CF the summation is on (u*1j)*r, in the second one its (phi -1j)1j*r. Maybe that's part of 'your' trick? $\endgroup$
    – Kermittfrog
    Apr 28, 2020 at 6:49
  • $\begingroup$ @Kermittfrog If you look at the integral in the first case, the argument in the characteristic function is u - 1j (look at the integral) and, in the characteristic function, the input gets multiplied by 1j, hence you have (u - 1j)*1j in the first case. In the second case, if you look at the integral, you should notice that phi=u. I just moved the -1j inside the other function because, well, I'm using a matrix as opposed to a function in the integral here. $\endgroup$
    – Stéphane
    Apr 28, 2020 at 15:08
  • $\begingroup$ @Kermittfrog The recursion is provided in the Heston and Nandi (2000) paper. It can also be found in the Christoffersen, Heston and Jacobs (2006) paper and, obviously, in his Toolbox as a MATLAB code. $\endgroup$
    – Stéphane
    Apr 28, 2020 at 15:11
  • $\begingroup$ Ok. Then, make sure that the char function is exact; and that your complex integration does not "bleed" any contributions from higher frequencies. If I remember correctly, that where the problems I used to have. But then, I switched over to FFT/FrFT after that... $\endgroup$
    – Kermittfrog
    Apr 29, 2020 at 7:00
  • $\begingroup$ I managed to reproduce the results of Christofference quite closely. It turns out that I needed a finer grid, especially close to zero (the integrand gets close to zero very fast, so finer points further away don't matter as much). For the moment, I have an ad hoc solution and it's going to be enough to finish my assignment. Thanks for the support. If you write up a summary as an answer, I'll be happy to pick it. $\endgroup$
    – Stéphane
    Apr 29, 2020 at 15:52

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.