Simulating GARCH Option Prices to Recover Black-Scholes Implied Volatility

# Black-Scholes implied volatility using a GARCH model


# Black-Scholes implied volatility using a GARCH model












Why I'm not getting the same Black-Scholes implied volatility values as the ones given in the paper "Asset pricing with second-order Esscher transforms" (2012) by Monfort and Pegoraro?

The authors consider the following model:

The stock return $$y_{t+1} = \mu_0 +\mu_1 y_t + \sigma_{t+1} \epsilon_{t+1} \quad \quad (1')$$

The variance process (GARCH) is given by $$\sigma^2_{t+1} = \omega_0 + \omega_1 (\sigma_t \epsilon_t - \omega_2)^2 + \omega_3 \sigma_t^2\quad \quad (1'')$$ where $$\epsilon_{t+1} - \epsilon_t \sim N(0,1) \quad \mbox{in} \quad \mathbb{P}\quad \quad (1''')$$

Under the risk-neutral probability the model becomes

The stock return

$$y_{t+1} = r_{t+1} -\frac{1}{2} (\sigma^\mathbb{Q}_{t+1})^2 + \sigma^\mathbb{Q}_{t+1} \epsilon^\mathbb{Q}_{t+1}\quad \quad (2')$$

The variance process (GARCH) is given by $$(\sigma^\mathbb{Q}_{t+1})^2 =\sigma^2_{t+1} \varphi \quad \phi:=e^{-\beta_0},\quad \beta_0\in \mathbb{R}\quad \quad (2'')$$

where

$$\epsilon^\mathbb{Q}_{t+1} - \epsilon^\mathbb{Q}_t \sim N(0,1) \quad \mbox{in} \quad \mathbb{Q}\quad \quad (2''')$$

Then the price of a European call option at date $t$ with underlying stock price $S_t$, residual maturity $\tau$ and moneyness strike $\kappa_t = \frac{K}{S_t}$ is given by \begin{align*} C(t,\tau,\kappa_t) &=\exp\bigg(-\sum_{i=1}^{\tau} r_{t+i} \bigg) E^\mathbb{Q}[S_{t+\tau}-K]^+\\ &\approx \exp\bigg(-\sum_{i=1}^{\tau} r_{t+i}\bigg) S_t \frac{1}{N} \sum_{n=1}^{N} \bigg[\exp(y^{(N)}_{t+1}+y^{(N)}_{t+2}+ \cdots + y^{(N)}_{t+\tau}) - \kappa_t\bigg]^+ \end{align*} where for each $n \in \{1,2,\cdots N\}$, the values $y^{(N)}_{t+1}+y^{(N)}_{t+2}+ \cdots + y^{(N)}_{t+\tau}$ entering in the pricing formula are simulated from the no arbitrage risk-neutral conditional distribution (2).

The authors provide the following

$\mu_0 = 0.00022$

$\mu_1 = 0.13899$

$\omega_0=4.54e-07$

$\omega_1 = 0.07105$

$\omega_2 = 0.00405$

$\omega_3 = 0.92353$

$r=\frac{0.05}{365}=0.000137$ daily interest rate

$N=10000$ number of simulations

$S_t=1$ initial stock price

$\kappa_t \in (0.8,1.2)$ moneyness

$\tau = 30$ days

$\beta_0 \in \{0, -0.25, -0.5, -0.75\}$

My goal

Using the above information I'm trying to reproduce the Black-Scholes implied volatility.

Take the simplest case where $\beta_0 = 0$. The authors implied volatility is the line in blue in the following graph

My attempt to reproduce the implied volatility using R

```
library(tseries)
library(tidyverse)
library(tidyquant)
library(fOptions)
library(timeSeries)

# variance estimated parameters
omega0 <- 4.54e-07
omega1 <- 0.07105
omega2 <- 0.00405
omega3 <- 0.92353
    
# long-run variance
sigma_squared_0 <- omega0/(1-omega1 * (1+omega2^2) - omega3)

var_par <- list(omega0=omega0, omega1=omega1, omega2=omega2, 
omega3=omega3, initial_var=initial_var)
    
# simulation and option info
r <- 0.05 # annualized
S0 <- 1
K <- seq(.8,1.2,.05)
days2maturity <- 30
    
sim <- 10000
    
option_par <- list(r=r, S0=S0, K=K, days2maturity=days2maturity, sim=sim)
    
beta0 <- c(0, -0.25, -0.5, -0.75)
    
set.seed(100)
    
main_cal <- function(var_par=var_par, option_par=option_par, beta=beta0[1]){
    t2m <- option_par$days2maturity
N <- option_par$days2maturity
    dt <- t2m / N # time step = 1
    time_stamp <- seq(0, t2m, dt)
    r <- option_par$r / 365 # interest rate per day
        
    sim <- option_par$sim
        
    # generate z (normal random variables doesn't depend on time)
    z <- matrix(stats::rnorm(sim*N, 0, 1), nrow = N)
        
    epsilon <- matrix(nrow=N+1, ncol = sim)
        
    epsilon[1,] <- 0
        
    sigma_squared <- matrix(nrow=N+1, ncol = sim)
        
    sigma_squared[1,] <- var_par$initial_var
        
    omega0 <- var_par$omega0
omega1 <- var_par$omega1
    omega2 <- var_par$omega2
omega3 <- var_par$omega3
        
    for (i in 1:N){
     z[i,] <- (z[i,] - mean(z[i,])) / sd(z[i,]) # ensure std norm
     epsilon[i+1,] <- epsilon[i,] + sqrt(dt) * z[i,]
        
     # simulate the variance process    
     sigma_squared[i+1,] <- omega0 
     + omega1*(sqrt(sigma_squared[i,]) * epsilon[i+1,] - omega2)^2 
     + omega3 * sigma_squared[i,]
    }
        
    phi <- exp(-beta) # = 1 when beta = 0
        
    sigma_squared_Q <- phi * sigma_squared
        
    # generate z^Q
    zQ <- matrix(stats::rnorm(sim*N, 0, phi), nrow = N)
        
    epsilonQ <- matrix(nrow=N+1, ncol = sim)
        
    epsilonQ[1,] <- 0
        
    S <- matrix(nrow=N+1, ncol = sim)
        
    S[1,] <- S0
        
    for (i in 1:N){
            
     zQ[i,] <- (zQ[i,] - mean(zQ[i,])) / sd(zQ[i,])
            
     epsilonQ[i+1,] <- epsilonQ[i,] + sqrt(dt) * zQ[i,]
            
     # stock formula
     S[i+1,] <- S[i,]*exp((r-0.5*sigma_squared_Q[i,])*dt 
     + sqrt(sigma_squared_Q[i,])*(epsilonQ[i+1,]-epsilonQ[i,]))

    }
        
    ST <- S[N+1,]
        
    optPrice <- c()
        
    for(k in 1:length(option_par$K)){
     payoff <- pmax(ST - K[k], 0)
            
     optPrice[k] <- exp(-r * t2m) * mean(payoff) 
    }
        
    result <- list(optPrice=optPrice,
    epsilon=epsilon,
    sigma_squared=sigma_squared,
    sigma_squared_Q=sigma_squared_Q,
    epsilonQ=epsilonQ,
    S=S)
        
    return(result)
        
}
    
data <- main_cal(var_par=var_par, option_par=option_par, beta=beta0[1])
    
optionPrice <- data$optPrice
    
optionPrice
    
sig <- data$sigma_squared
    
S <- data$S
    
imp_vol <- c()
    
for(i in 1:length(K)){
 fun <- function(x){
    fOptions::GBSOption(TypeFlag = "c", S = S0, X = K[i], 
    Time = days2maturity/365, r =r, sigma = x, b = r)@price 
    - optPrice[i]
    }
        
    imp_vol[i] <- stats::uniroot(fun, c(1e-04, 2), tol =1e-06)$root
}
    
plot(K, imp_vol, type="l")
```

However, I'm having this error

```
Error in stats::uniroot(fun, c(1e-04, 2), tol = 1e-06) : 
f() values at end points not of opposite sign
```

Also, investigating the plot, I found that to reproduce it using the above algorithm the option prices must be around the following values

```
plot_imp_vol <- c(.26,.241,.228,.222,.219,.223,.225,.23,.237)
        
optPrice <- c(0.2048832, 0.1556492, 0.1079794, 0.06544197, 
0.03310882, 0.01421061, 0.00502847, 0.001603875, 0.0004927986)
```

which are different than mine.

I would appreciate the help in pointing out what I'm missing and what should I adjust to obtain the same plot.

Thank you, paper