I am trying to replicate Duan's results from his 1995 Paper, "The GARCH Option Pricing Model". I have written this code in Python myself, and using his parameters I consistently seem to obtain results significantly below his results. As an example, if I run the code with 30 days as Time to Maturity of the Option and number of simulations being 50,000, I obtain a result of 260.162, whereas Duan reports results of 266.75. I was wondering if anybody who was familiar with GARCH Option Pricing could assist me in my troubles. Thankyou in advance. (Note that results are reported as 10,000 times, if anyone was wondering how a call option with initial stock price of 1, strike price of 1, t = 30 days and r = 0 could be price as 266.75).
import numpy as np beta0 = 0.00001524 #GARCH Parameter Omega beta1 = 0.7162 #GARCH Parameter associated with lagged variance term beta2 = 0.1883 #GARCH Parameter associated with lagged innovation lamba = 0.007452 #RiskNeutral Parameter H0 = beta0/(1-(1+lamba)*beta2-beta1) #initial variance #H0 = (0.8)**2*(omega/(1-(1+lamba)*alpha-beta)) #initial conditional variance #H0 = (1.2)**2*(omega/(1-(1+lamba)*alpha-beta)) #initial conditional variance S0 = 1 #initial stock price K = 1 #strike price r = 0 #risk-free interest rate t = 0 #start time Td = 30 #time in days i = 50000 #number of simulations discount_factor = np.exp(-r*(Td/365)) #discount factor dt = np.dtype(np.float16) h = np.zeros([i,Td], dtype = dt) e = np.zeros([i,Td], dtype = dt) t = range(0,Td,1) S = S0*np.ones([i], dtype = dt) DH = np.ones([i], dtype = dt) z = np.random.standard_normal([i,Td]) for x in range(0,i-1): h[x,0] = H0 e[x,0] = H0*np.random.normal(0,1,) for y in range(0,Td-1): h[x,y+1] = beta0 + h[x,y]*(beta1 + beta2*(z[x,y] - lamba)**2) e[x,y+1] = np.sqrt(h[x,y+1])*z[x,y+1] sumh = h.sum(axis=1) sume = e.sum(axis=1) for x in range(0,i): S[x] = S[x]*np.exp((Td)*(r/365) - 0.5*sumh[x] + sume[x]) DH[x] = np.exp(-0.5*sumh[x] + sume[x]) for x in range(0,i): if S[x]>=K: DH[x] = DH[x] else: DH[x] = 0 for x in range(0,i): S[x] = np.maximum(S[x] - K, 0) S[x] = discount_factor*S[x] CallPrice = np.average(S) StdDev = np.std(S) UpperLimit = CallPrice + 1.96*StdDev/i LowerLimit = CallPrice - 1.96*StdDev/i DeltaHedge = np.average(DH) #print("The Call Price is:", CallPrice*10000) print("The Call Price is:", CallPrice*10000) print("The lower bound of the 95% confidence interval is:", LowerLimit) print("The upper bound of the 95% confidence interval is:", UpperLimit) print("The corresponding Option Delta is:", DeltaHedge)