Reproducing Black–Scholes and GARCH Option Prices in Duan’s Model

# Problem matching prices of Black-Scholes vs. GARCH(1,1) in Duan (1995)


# Problem matching prices of Black-Scholes vs. GARCH(1,1) in Duan (1995)












In the paper of Duan (1995) the author compare European call option prices using Black-Scholes model vs. GARCH(1,1)-M model (GARCH-in-mean). To be brief, the author fits the following GARCH(1,1)-M model

$$\ln(X_t/X_{t-1}) = r + \lambda \sqrt{h_t} - \frac{1}{2} h_t + \varepsilon_t $$ $$ \varepsilon_t |\phi_{t-1} \sim N(0,h_t) \quad \mbox{under measure P} $$ $$ h_t = \alpha_0 + \alpha_1 \varepsilon_{t-1}^2 + \beta_1 h_{t-1} \tag{1}$$

to S&P 100 daily index series from `1986-01-02` to `1989-12-15`, where $r$ is the risk-free interest rate, $\lambda$ is the constant unit risk premium and $\phi_t$ is the information set. For GARCH(1,1), $X_t$ and $h_t$ together serve as the sufficient statistics for $\phi_t$.

The author also provides the risk-neutral dynamics of (1) as follows $$\ln(X_t/X_{t-1}) = r - \frac{1}{2} h_t + \xi_t $$ $$ \xi_t |\phi_{t-1} \sim N(0,h_t) \quad \mbox{under measure Q} $$ $$ h_t = \alpha_0 + \alpha_1 (\xi_{t-1}-\lambda \sqrt{h_{t-1}})^2 + \beta_1 h_{t-1} \tag{2}$$

They also state that the price of a European call option with exercise price $K$ is given by $$C^{GH}_0 = e^{-rT} \mathbb{E}^Q \bigg[ \max(X_T - K, 0) \bigg]$$ where the terminal asset price is derived as $$X_T = X_0 \exp \bigg[rT -\frac{1}{2} \sum_{s=1}^T h_s + \sum_{s=1}^T \xi_s \bigg] \tag{3} $$

For the computation part, the author provides the estimated GARCH parameters and some other assumptions used in their calculations. Their assumptions and results are in the following R code

```
#----------------------------------------------------
# information and results given by the author
#----------------------------------------------------
# estimated parameters of the GARCH(1,1)-M model
alpha0 <- 1.524e-05
alpha1 <- 0.1883
beta1 <- 0.7162
lambda <- 7.452e-03

r <- 0 # interest rate
K <- 1 # strike price is set at $1 (stated like this in the paper)
T2M <- c(30, 90, 180)  # days to maturity
S_over_K <-c(.8, .9, .95, 1, 1.05, 1.1, 1.2)  # moneyness ratios considered (stated like this in the paper)
Sigma <- 0.2413 / sqrt(365)

simulations <- 5*10^4  # Monte Carlo simulations

maturity <- c(rep(T2M[1], length(S_over_K)), rep(T2M[2], length(S_over_K)), rep(T2M[3], length(S_over_K)))
SK_ratio <- rep(S_over_K, length(T2M))

# Black-Scholes option prices
BS_price <- c(.1027, 18.238,89.79, 276.11, 600.41, 1027.9, 2001, 13.016, 118.54, 257.79, 478, 779.68, 1152.2, 2036.6, 66.446, 261.3, 438.17, 675.54, 970.49, 1318, 2133.4)

# GARCH prices when sqrt(h_1) = Sigma
GH_price <- c(.9495, 20.93, 86.028, 266.75, 596.13, 1030.2, 2003.5, 15.759, 116.06, 251.02, 468.9, 772.35, 1149.4, 2040.5, 68.357, 257.09, 431.7, 668.5, 964.29, 1313.9, 2134.7)

df <- data.frame(maturity, SK_ratio, BS_price, GH_price)
df
```

I'm trying to get their results by using the following code.

```
#-------------------------------------------------------
# my attempt
#-------------------------------------------------------
library(rugarch)
library(stats)
library(quantmod)
library(PerformanceAnalytics)
library(fOptions)
# download S&P100 daily index series
fromDate <- "1986-01-02"
toDate <- "1989-12-15"
quantmod::getSymbols("^OEX", from = fromDate, to = toDate)

# calculate the log returns
sp100ret <- na.omit(PerformanceAnalytics::CalculateReturns(OEX$OEX.Close, method = "log"))

# set GARCH(1,1)-M model specification
garchspec <- ugarchspec(mean.model = list(armaOrder = c(0, 0), include.mean=TRUE, archm=TRUE, archpow=1),
                        variance.model = list(model = "sGARCH", garchOrder=c(1,1)), 
                        distribution.model = "norm")
# fit the model to the daily returns data
sp100Fit <- ugarchfit(data = sp100ret, spec = garchspec)

# output the GARCH estimated parameters
garchcoef <- coef(sp100Fit)
garchcoef

# assume the initial conditional variance = stationary level i.e, sqrt(h_1) = Sigma
h_1 <- alpha0/(1-alpha1 - beta1)

S0 <- S_over_K * K

my_BS_price <- c()

# calculate Black-Scholes prices
for (i in 1:3){
  my_BS_price <- c(my_BS_price, fOptions::GBSOption(TypeFlag="c", S=S0, X=K, Time=T2M[i], r=r, sigma=h_1, b=r)@price)
}

# create empty matricies for the standard normal random variables under Q
epsilonQ_T30 <- matrix(stats::rnorm((T2M[1]+1)*simulations, 0, 1), nrow = T2M[1]+1, ncol= simulations)
epsilonQ_T90 <- matrix(stats::rnorm((T2M[2]+1)*simulations, 0, 1), nrow = T2M[2]+1, ncol= simulations)
epsilonQ_T180 <- matrix(stats::rnorm((T2M[3]+1)*simulations, 0, 1), nrow = T2M[3]+1, ncol= simulations)

epsilonQ_T30[1,] <- 0
epsilonQ_T90[1,] <- 0
epsilonQ_T180[1,] <- 0

h_T30 <- matrix(nrow = T2M[1]+1, ncol = simulations)
h_T90 <- matrix(nrow = T2M[2]+1, ncol = simulations)
h_T180 <- matrix(nrow = T2M[3]+1, ncol = simulations)

h_T30[1,] <- h_1
h_T90[1,] <- h_1
h_T180[1,] <- h_1

# for loop to calculate the conditional variance process Eq(2) maturity 30
for (t in 2:31){
  h_T30[t, ] <- alpha0 + alpha1 * h_T30[t-1,] * (epsilonQ_T30[t-1,] - lambda)^2 + beta1 * h_T30[t-1,]
  epsilonQ_T30[t, ] <- epsilonQ_T30[t, ] * sqrt(t) + lambda * t # risk-neutral
}

# for loop to calculate the conditional variance process Eq(2) maturity 90
for (t in 2:91){
  h_T90[t, ] <- alpha0 + alpha1 * h_T90[t-1,] * (epsilonQ_T90[t-1,] - lambda)^2 + beta1 * h_T90[t-1,]
  epsilonQ_T90[t, ] <- epsilonQ_T90[t, ] * sqrt(t) + lambda * t # risk-neutral
}

# for loop to calculate the conditional variance process Eq(2) maturity 180
for (t in 2:181){
  h_T180[t, ] <- alpha0 + alpha1 * h_T180[t-1,] * (epsilonQ_T180[t-1,] - lambda)^2 + beta1 * h_T180[t-1,]
  epsilonQ_T180[t, ] <- epsilonQ_T180[t, ] * sqrt(t) + lambda * t # risk-neutral
}

# remove first row
h_T30 <- h_T30[-c(1), ]
h_T90 <- h_T90[-c(1), ]
h_T180 <- h_T180[-c(1), ]

epsilonQ_T30 <- epsilonQ_T30[-c(1),]
epsilonQ_T90 <- epsilonQ_T90[-c(1),]
epsilonQ_T180 <- epsilonQ_T180[-c(1),]

# calculate the exponential part in terminal value of stock Eq(3)
ST_expPart_T30 <- exp(r*T2M[1] - 0.5 * colSums(h_T30) + colSums(sqrt(h_T30)) * colSums(epsilonQ_T30))
ST_expPart_T90 <- exp(r*T2M[2] - 0.5 * colSums(h_T90) + colSums(sqrt(h_T90)) * colSums(epsilonQ_T90))
ST_expPart_T180 <- exp(r*T2M[3] - 0.5 * colSums(h_T180) + colSums(sqrt(h_T180)) * colSums(epsilonQ_T180))

# calculate the terminal value of the stock for each initial value Eq(3)
ST_T30 <- matrix(nrow = length(S0), ncol=simulations)
ST_T90 <- matrix(nrow = length(S0), ncol=simulations)
ST_T180 <- matrix(nrow = length(S0), ncol=simulations)

for(i in 1:length(S0)){
  ST_T30[i,] <- S0[i] * ST_expPart_T30
  ST_T90[i,] <- S0[i] * ST_expPart_T90
  ST_T180[i,] <- S0[i] * ST_expPart_T180
}

# calculate expectation of the discounted call option payoff
C_T30 <- c()
C_T90 <- c()
C_T180 <- c()

for(i in 1:length(S0)){
  C_T30[i] <- exp(-r*T2M[1]) * mean( pmax(ST_T30[i,] -K, 0) )
  C_T90[i] <- exp(-r*T2M[2]) * mean( pmax(ST_T90[i,] -K, 0) )
  C_T180[i] <- exp(-r*T2M[3]) * mean( pmax(ST_T180[i,] -K, 0) )
}

my_GH_price <- c(C_T30, C_T90, C_T180)

my_df <- data.frame(maturity, SK_ratio, my_BS_price, my_GH_price)
my_df
```

where $$\xi_t = \sqrt{h_t} \varepsilon^Q_t ;\quad \varepsilon^Q_t \sim N(0,t) \quad \mbox{under measure Q or }\quad \varepsilon^Q_t \sim N(\lambda t,t) \quad \mbox{under measure P}$$

I have three problems.

- I didn't get the exact estimated parameters as theirs. Is the `ugarchspec` wrong?

- I used the author's estimated parameters to calculate option prices using Black-Scholes formula, but I didn't get the same numbers. I don't understand why the BS-price is not the same. Any idea?

- I didn't get the same numbers as the paper for the GARCH model prices either. Is there anything wrong with how I modeled the risk-neutral dynamics of (2)?

I appreciate the help.