Gauss-Hermite Quadrature

Efficient numerical integration method with weight function exp(-x^2)

You need this for implementing Kennedy's method

There are two versions:

Probabilists’ Gauss-Hermite module: integration weight is the standard normal PDF: exp(-x^2/2)
Physicists’ Gauss-Hermite module: integration weight is exp(-x^2)

We mostly use Probabilists’ Gauss-Hermite module. You still need to devide the weight by sqrt(2*pi) for use with normal PDF



In [1]:

    
import numpy as np
import numpy.polynomial as nppoly
import scipy
import scipy.stats as ss
import matplotlib.pyplot as plt



In [2]:

    
const = 1/np.sqrt(2.0*np.pi)



In [4]:

    
z, w = nppoly.hermite_e.hermegauss(deg=20)
w = w*const
print(z)
print(w)









    



[-7.61904854 -6.51059016 -5.57873881 -4.73458133 -3.94396735 -3.18901482
 -2.45866361 -1.74524732 -1.04294535 -0.34696416  0.34696416  1.04294535
  1.74524732  2.45866361  3.18901482  3.94396735  4.73458133  5.57873881
  6.51059016  7.61904854]
[  1.25780067e-13   2.48206236e-10   6.12749026e-08   4.40212109e-06
   1.28826280e-04   1.83010313e-03   1.39978374e-02   6.15063721e-02
   1.61739334e-01   2.60793063e-01   2.60793063e-01   1.61739334e-01
   6.15063721e-02   1.39978374e-02   1.83010313e-03   1.28826280e-04
   4.40212109e-06   6.12749026e-08   2.48206236e-10   1.25780067e-13]



In [5]:

    
pdf = ss.norm.pdf(z)
plt.plot(z, np.log(w))
plt.plot(z, np.log(pdf))
plt.grid()
plt.show()

Exact integration of polynomials with degree upto 2deg-1*



In [14]:

    
z, w = nppoly.hermite_e.hermegauss(deg=3)
w = w*const
sum(w)









    Out[14]:





1.0



In [15]:

    
# Let's test on the moments of normal distribution
deg = np.array([2,4,6,8,10,12,14])
moments = [sum(z**2 * w), sum(z**4 * w), sum(z**6 * w), sum(z**8 * w), sum(z**10 * w), sum(z**12 * w), sum(z**14 * w)]
print(moments)









    



[0.99999999999999978, 2.9999999999999996, 9.0, 27.000000000000004, 81.000000000000028, 243.00000000000009, 729.00000000000045]



In [16]:

    
# luckily we know the exact answer: (2*deg-1)!! 
scipy.misc.factorial2([1,3,5,7,9,11,13])









    Out[16]:





array([  1.00000000e+00,   3.00000000e+00,   1.50000000e+01,
         1.05000000e+02,   9.45000000e+02,   1.03950000e+04,
         1.35135000e+05])



In [17]:

    
# Find out upto which degree integration is correct
deg[np.abs(moments - scipy.misc.factorial2([1,3,5,7,9,11,13])) < 0.1 ]









    Out[17]:





array([2, 4])

Overall GHQ is very accurate for integrating smooth functions

Let's test on Geometric Brownian Motion:

$ S_T = S_0 exp\left(\sigma\sqrt{T} z - \frac12 \sigma^2 T\right)$



In [18]:

    
spot = 100
texp = 2
vol = 0.2



In [19]:

    
z = np.linspace(-5,5,10)
price = spot * np.exp(vol*np.sqrt(texp)*z - 0.5*vol*vol*texp)
print(price)









    



[  23.35839909   31.983697     43.79396334   59.96527623   82.10799112
  112.42710165  153.94181507  210.78620794  288.6209016   395.19675245]



In [ ]:



In [20]:

    
# Let's check the expectation of the prices are same as 100 (assuming 0 interest rate)
z, w = nppoly.hermite_e.hermegauss(deg=10)
w = w*const
price = spot * np.exp(vol*np.sqrt(texp)*z - 0.5*vol*vol*texp)
price_mean = sum(price * w)
price_mean - 100









    Out[20]:





-1.4210854715202004e-14



In [21]:

    
plt.plot(price, w, 'o-')
plt.grid()
plt.show()



In [ ]:

Gauss-Hermite Quadrature

Efficient numerical integration method with weight function exp(-x^2)

You need this for implementing Kennedy's method

Exact integration of polynomials with degree upto 2*deg-1

Overall GHQ is very accurate for integrating smooth functions

Exact integration of polynomials with degree upto 2deg-1*