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,