Replicating Duan’s GARCH Option Pricing with Monte Carlo Simulation

# GARCH Option Pricing Model (Duan 1995)


# GARCH Option Pricing Model (Duan 1995)












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