Volatility Smile and Stochastic Volatility Model

i. Introduction

Option is a financial derivative that its value depends on the other asset. There are a few methods that are widely adopted in option pricing:

• Lattice Models(Binomial Tree, Trinomial Tree,…)
• Analytical Models(Black-Scholes Model,…)
• Monte-Carlo Simulation Methods
• Finite Difference Methods

In this project, analytical pricing method(Black-Scholes) and the arousing issues are discussed, and Stochastic Volatility Model is introduced as a better way to fix the arousing problem.

ii.Black-Scholes Pricing Model

By Black-Scholes Model, European options can be calculated by

\(C= Se^{-q(T-t)}N(d_{1})-K e^{-R(T-t)} N(d_{2})\),

\(P= K e^{-R(T-t)} N(-d_{2})- Se^{-q(T-t)}N(-d_{1}), where\)

\(d_{1}= \frac{log(\frac{S}{K})+(R-q+\frac{\sigma^{2}}{2})(T-t)}{\sigma \sqrt{T-t}},d_{2}= d_{1}- \sigma \sqrt{T-t}\)

C is call option price, P is put option price, S is the underlying stock price, K is the strike price, R is the risk-free rate, T-t is time to maturity, q is the dividend rate and N(.) is cumulative standard normal distribution.

There are a few assumptions under Black-Scholes Model:

1. Risk-free(R) interest rate is constant
2. Drift rate(\(\mu\)) is constant.
3. Volatility(\(\sigma\)) is constant.
4. The underlying asset price(S) follows a geometric Brownian Motion.

For simplicity, we only focus on the plain vanilla type European option(without dividend).

For a specific European option, almost every parameters can be directly observed by the market, except for only one parameter- volatility. We can make use of those observable parameters and Black-Scholes Model to calculate volatility, which is frequently referred to "Implied Volatility".

iii. Volatility Skew and Smile

The closed form solution from Black and Scholes provides a quick and easy way to compute the price of European options. However, the assumption of constant volatility arouses a issue: volatility skew(or volatility smile).

callbs <- function(S, K, R, T, sig, q=0) {
  d1 <- (log(S/K) + (R +q+ sig^2/2)*T) / (sig*sqrt(T))
  d2 <- d1 - sig*sqrt(T)
  c <- S*exp(-q*T)*pnorm(d1) - K*exp(-R*T)*pnorm(d2)
  if(c>=0){return(c)}
  else{return(0)}
}
put_bs <- function(S, K, R, T, sig, q=0) {
  d1 <- (log(S/K) + (R + q+ sig^2/2)*T) / (sig*sqrt(T))
  d2 <- d1 - sig*sqrt(T)
  p <- K*exp(-R*T)*pnorm(-d2)-S*exp(-q*T)*pnorm(-d1)
  if(p>=0){return(p)}
  else{return(0)}
}

suppressPackageStartupMessages(library(quantmod))
suppressMessages(suppressWarnings(getSymbols("AAPL",from= Sys.Date()-3, auto.assign = TRUE)))
S= as.numeric(AAPL$AAPL.Close)[1]

suppressPackageStartupMessages(library(quantmod))
AAPL.OPTS <- getOptionChain("AAPL",NULL)
AAPL.OPTS <- AAPL.OPTS[[length(AAPL.OPTS)]]


#call
call <- AAPL.OPTS$calls
C= call$Last
K= call$Strike
R= 0.0013
# US treasury bill 2 year rate= 0.13%


year= paste0("20",substring(unlist(labels(call))[1],5,6))
month= substring(unlist(labels(call))[1],7,8)
day= substring(unlist(labels(call))[1],9,10)
date= paste0(year,"-",month,"-",day)
T= as.numeric(difftime(date,Sys.Date(), units = "days")/365)
  
sig_impl <- function(S, K, R, T, C){
  High= 1
  Low=0
  while((High-Low)>0.000001){
    if(callbs(S,K,R,T,(High+Low)/2) > C){
      High= (High+Low)/2
    }else{
      Low= (High+Low)/2
    }
  }
  return((High + Low)/2)}

Impl_Vol=c()
for(i in 1:nrow(call)){
  Impl_Vol[i]= sig_impl(S, K=call$Strike[i],R,T,C=call$Last[i])
}
#print(Impl_Vol)
#plot(x= call$Strike, y=Impl_Vol)

df_c <- as.data.frame(cbind(call$Strike, Impl_Vol))
names(df_c)[1] <- "Strike"
suppressPackageStartupMessages(library(ggplot2))
g1 <- ggplot(data= df_c)+ geom_point(aes(x= Strike, y= Impl_Vol))+ xlab("Strike Price")+ ylab("Implied Volatility")+ ggtitle("Volatility Skew for Call options")

suppressPackageStartupMessages(library(quantmod))
AAPL.OPTS <- getOptionChain("AAPL",NULL)
AAPL.OPTS <- AAPL.OPTS[[length(AAPL.OPTS)]]

#puts
put <- AAPL.OPTS$puts
P= put$Last
K= put$Strike
R= 0.0013
# US treasury bill 2 year rate= 0.13%


year= paste0("20",substring(unlist(labels(put))[1],5,6))
month= substring(unlist(labels(put))[1],7,8)
day= substring(unlist(labels(put))[1],9,10)
date= paste0(year,"-",month,"-",day)
T= as.numeric(difftime(date,Sys.Date(), units = "days")/365)
  
sig_impl <- function(S, K, R, T, P){
  High= 1
  Low=0
  while((High-Low)>0.000001){
    if(put_bs(S,K,R,T,(High+Low)/2) > P){
      High= (High+Low)/2
    }else{
      Low= (High+Low)/2
    }
  }
  return((High + Low)/2)}

Impl_Vol=c()
for(i in 1:nrow(put)){
  Impl_Vol[i]= sig_impl(S, K=put$Strike[i],R,T,P=put$Last[i])
}
#print(Impl_Vol)
#plot(x= call$Strike, y=Impl_Vol)

df_p <- as.data.frame(cbind(put$Strike, Impl_Vol))
names(df_p)[1] <- "Strike"
suppressPackageStartupMessages(library(ggplot2))

g2 <- ggplot(data= df_p)+ geom_point(aes(x= Strike, y= Impl_Vol))+ xlab("Strike Price")+ ylab("Implied Volatility")+ ggtitle("Volatility Skew for Put options")

library(gridExtra)
grid.arrange(g1,g2,ncol=2)

iv. Stochastic Volatility Model

Stochastic Volatility Model explains the reason of volatility skew/smile. Stochastic Volatility Model assumes the volatility on day t, is partially determined by some unpredicted factors.

\(dS_{t}=\mu S_{t} dt+ \sqrt {v_{t}} S_{t} dZ_{1}\)

\(dv_{t}= \alpha(S_{t},v_{t},t)dt+\sigma\beta(S_{t}, v_{t}, t) \sqrt v_{t} dZ_{2}\)

\(dZ_{1} \cdot dZ_{2}= \rho dt\)

, where \(\mu\) is the drift rate of the stock returns, \(\sigma\) is the volatility of volatility, \(\rho\) is the correlation between the random change of stock returns and the random change in volatility, and \(dZ_{1}\) and \(dZ_{2}\) are Brownian motions.

Consider forming a portfolio \(\Pi\) by shorting \(\Delta\) shares of stocks with initial stock price S, shorting \(\Delta_{1}\) units of another asset whose value \(V_{1}\) depends on volatility, and long an option with value of \(V(S,v,t)\). Then the portfolio value \(\Pi= V-\Delta S- \Delta_{1} V_{1}\),

\(d\Pi=\{\frac{\partial V}{\partial t}+\frac {1}{2} vS^{2} \frac{\partial ^{2}V}{\partial S^{2}}+\rho \sigma v \beta S \frac {\partial^{2} V}{\partial v \partial S}+\frac{1} {2} \sigma^{2}v \beta^{2} \frac{\partial^{2} V}{\partial v^{2}}\}dt-\Delta_{1} \{\frac{\partial V_{1}}{\partial t}+\frac{1}{2}vS^{2} \frac{\partial^{2}V_{1}}{\partial S^{2}}+\rho \sigma v \beta S \frac{\partial^{2}V_{1}}{\partial v \partial S}+\frac{1}{2} \sigma^{2} v \beta^{2} \frac{\partial^{2} V_{1}}{\partial v^{2}}\}dt+\)

\(\{\frac{\partial V}{\partial S}-\Delta_{1} \frac{\partial V_{1}}{\partial S}-\Delta \}dS+\{\frac{\partial V}{\partial v}-\Delta_{1} \frac{\partial V_{1}}{\partial v}\}dv\)

Under no-arbitrage assumption, in order to make sure the portfolio \(\Pi\) only earns risk-free interest rate, the following equations should be satisfied:

\(\frac{\partial V}{\partial S}-\Delta_{1} \frac{\partial V_{1}}{\partial S} - \Delta =0\)

\(\frac {\partial V}{\partial v} -\Delta_{1} \frac{\partial V_{1}}{\partial v}=0\)

Thus the portfolio value will not depend on the change in stock price and the change in volatility. Then

\(d \Pi= \{\frac{\partial V}{\partial t}+\frac{1}{2} vS^{2} \frac{\partial^{2} V}{\partial S^{2}}+ \rho\sigma v \beta S \frac{\partial^{2} V}{\partial v \partial S}+\frac{1}{2} \sigma^{2} v \beta^{2} \frac{\partial^{2} V}{\partial v^{2}}\}dt-\Delta_{1}\{\frac{\partial V_{1}}{\partial t}+\frac{1}{2} vS^{2} \frac{\partial^{2} V_{1}}{\partial S^{2}}+\rho \sigma v \beta S \frac{\partial^{2} V_{1}}{\partial v \partial s}+\frac{1}{2} \sigma^{2} v \beta^{2} \frac{\partial^{2} V_{1}}{\partial v^{2}}\}dt\)

\(= r \Pi dt= r(V-\Delta S- \Delta_{1} V_{1})dt\)

Put all \(V\) terms on left hand side, and all \(V_{1}\) terms on right hand side, we get:

\(\frac{\frac{\partial V}{\partial t}+\frac{1}{2} v S^{2} \frac{\partial^{2} V}{\partial S^{2}}+\rho \sigma v \beta S \frac{\partial^{2}V}{\partial v \partial S}+\frac{1}{2} \sigma^{2} v \beta^{2} \frac{\partial^{2} V}{\partial v^{2}}+rS \frac{\partial V}{\partial S}-rV}{\frac{\partial V}{\partial v}} = \frac{\frac{\partial V_{1}}{\partial t}+ \frac{1}{2}v S^{2} \frac{\partial^{2} V_{1}}{\partial S^{2}}+ \rho \sigma v \beta S \frac{\partial^{2} V_{1}}{\partial v \partial S}+ \frac{1}{2} \sigma^{2} v \beta^{2} \frac{\partial^{2} V_{1}}{\partial v^{2}}+ rS \frac{\partial V_{1}}{\partial S} - rV_{1}}{\frac{\partial V_{1}}{\partial v}}\)

There must be some function \(f(S,v,t)\) that satisfies the equation above.

\(\frac{\partial V}{\partial t}+ \frac{1}{2}v S^{2} \frac{\partial^{2} V}{\partial S^{2}}+ \rho \sigma v \beta S \frac{\partial^{2} V}{\partial v \partial S}+ \frac{1}{2} \sigma^{2} v \beta^{2} \frac{\partial^{2} V}{\partial v^{2}}+ rS \frac{\partial V}{\partial S}- rV = f(S,v,t) \frac{\partial V}{\partial v}=-(\alpha - \phi \beta \sqrt v) \frac{\partial V}{\partial v}\)

where the function \(f(S,t,v)\) has been set as \(-(\alpha- \phi \beta \sqrt v) \frac{\partial V}{\partial v}\), \(\alpha\) is drift and \(\beta\) is diffusivity.

v. European Option Pricing by The Heston Model

The Heston model pricing fromula is consistent with Black-Scholes Model, we can write the call option pricing formula as

\(C(S,v,t)= SP_{1}-Ke^{-R(T-\tau)}P_{2}\)

According to Steven Heston(1993), the solution of the pricing method should be:

Define \(\tau=T-t\) is defined as time to maturity, \(C_{i}= \frac{\partial C}{\partial i}, i=x, v\), and \(C_{ij}= \frac{\partial^{2} C}{\partial i \partial j}, i,j=x,v.\), let \(x=ln(S)\), subsitute \(C\) into the SDE derived from The Heston Model, then \(P_{0}\) and \(P_{1}\) must satify:

\(-\frac{\partial C}{\partial \tau}+\frac{1}{2}vC_{xx}-\frac{1}{2}vC_{x}+\frac{1}{2}\sigma^{2}vC_{vv}+\rho \sigma v C_{xv}-\kappa (v-\theta) C_{v}=0\)

\(=> \frac{1}{2}v\frac{\partial^{2}P_{j}}{\partial x^{2}}+\rho \sigma v \frac{\partial^{2}P_{j}}{\partial x \partial v}+\frac{1}{2} \sigma^{2} v \frac{\partial^{2} P_{j}}{\partial v^{2}}+ (r+u_{j}v)\frac{\partial P_{j}}{\partial x}+(a_{j}-b_{j}v)\frac{\partial P_{j}}{\partial v}+\frac{\partial P_{j}}{\partial t}=0, j=1,2\), where

\(u_{1}=\frac{1}{2}, u_{2}=-\frac{1}{2}, a= \kappa \theta, b_{1}=\kappa + \lambda- \rho \sigma, b_{2}= \kappa+ \lambda\)

Solving the differential equations by Fourier Transformation,

\(P_{j}=\frac{1}{2}+\frac{1}{\pi}\int_{0}^{\infty} Re\{\frac{e^{-i \phi lnK} f_{j}(x,v,\tau;\phi)}{i \phi}\}d\phi, j=1,2\), where

\(f_{j}(x,v,\tau;\phi)=e^{C_{j}(\tau;\phi)+D_{j}(\tau;\phi)v+i \phi x}\)

\(C(\tau;\phi)= r \phi i\tau+ \frac{a}{\sigma^{2}}\{(b_{j}-\rho \sigma \phi i+d)\tau-2ln(\frac{1- ge^{d\tau}}{1-g})\}\)

\(D(\tau;\phi)= \frac{b_{j}-\rho \sigma \phi i+d}{\sigma^{2}}(\frac{1-e^{d\tau}}{1-ge^{d\tau}})\)

\(g=\frac{b_{j}-\rho\sigma \phi i +d}{b_{j}-\rho \sigma \phi i-d}\)

\(d= \sqrt{(\rho \sigma \phi i- b_{j})^{2}- \sigma^{2}(2u_{j}\phi i-\phi^{2})}\)

vi. Results from Black-Scholes and Heston Model

We take the 12.May.2020 market quotes and data of APPLE.Inc(AAPL) from Yahoo Finance. Call option price on APPLE.Inc(AAPL) can be calculated by Black-Scholes Model and Heston SDE.

AAPL2 <- getSymbols("AAPL",auto.assign = FALSE, from= Sys.Date()-180)
ret <- diff(log(AAPL2$AAPL.Close))
ret[1] <- 0
vol <- sd(ret)
callbs(311.41, 120, 0.0013, 2.095776, vol, 0.0106)

AAPL2 <- getSymbols("AAPL",auto.assign = FALSE, from= Sys.Date()-180)
ret <- diff(log(AAPL2$AAPL.Close))
ret[1] <- 0
vol <- sd(ret)
callbs(311.41, 485, 0.0013, 2.095776, vol, 0.0106)

target <- 192.40
v0 = 0.2
vT = 0.2
rho= 0.5
k= 0.2
sig= 0.05

library(NMOF)
x <- c(v0,vT,rho,k,sig)
optim_func <- function(x){
  v0= x[1]
  vT= x[2]
  rho= x[3]
  k= x[4]
  sig= x[5]
  mse <- (callHestoncf(315.01,120,2.095776, 0.0013,0.0106, v0, vT, rho,k,sig)-target)^2
  return(mse)
}

optim(x, optim_func)

mat <- matrix(nrow=5, ncol=3)
mat[1,] <- c("v0","current variance",0.20940146)
mat[2,] <- c("vT", "long-run variance", 0.21366057)
mat[3,] <- c("rho","correlation between spot and variance", 0.50481539)
mat[4,] <- c("k","speed of mean-reversion", 0.21543664)
mat[5,] <- c("sigma", "volatility of variance", 0.04229108)
df <- as.data.frame(mat)
colnames(df) <- c("Parameters","Meaning","Value")
print(df)

##   Parameters                               Meaning      Value
## 1         v0                      current variance 0.20940146
## 2         vT                     long-run variance 0.21366057
## 3        rho correlation between spot and variance 0.50481539
## 4          k               speed of mean-reversion 0.21543664
## 5      sigma                volatility of variance 0.04229108

callHestoncf(311.41 ,120,2.095776, 0.0013,0.0106,0.20940146, 0.21366057, 0.50481539, 0.21543664, 0.04229108)

target <- 12.36
v0 = 0.2
vT = 0.2
rho= 0.5
k= 0.2
sig= 0.05

library(NMOF)
x <- c(v0,vT,rho,k,sig)
optim_func <- function(x){
  v0= x[1]
  vT= x[2]
  rho= x[3]
  k= x[4]
  sig= x[5]
  mse <- (callHestoncf(315.01,480,2.095776, 0.0013,0.0106, v0, vT, rho,k,sig)-target)^2
  return(mse)
}

optim(x, optim_func)

mat <- matrix(nrow=5, ncol=3)
mat[1,] <- c("v0","current variance",0.03401212)
mat[2,] <- c("vT", "long-run variance", 0.19923177)
mat[3,] <- c("rho","correlation between spot and variance", 0.54979724)
mat[4,] <- c("k","speed of mean-reversion", 0.30583280)
mat[5,] <- c("sigma", "volatility of variance", 0.08600963)
df <- as.data.frame(mat)
colnames(df) <- c("Parameters","Meaning","Value")
print(df)

##   Parameters                               Meaning      Value
## 1         v0                      current variance 0.03401212
## 2         vT                     long-run variance 0.19923177
## 3        rho correlation between spot and variance 0.54979724
## 4          k               speed of mean-reversion  0.3058328
## 5      sigma                volatility of variance 0.08600963

callHestoncf(311.41 ,485,2.095776, 0.0013,0.0106,0.03401212, 0.19923177, 0.54979724, 0.30583280, 0.08600963)

vii. Comparison and Conclusion

The actual call option market quotes and the calculated results are summarized as below:

mat <- matrix(nrow=2, ncol=4)
mat[1,] <- c("In-the-Money",184.8947,189.0168, 192.4)
mat[2,] <- c("Out-of-Money",0,11.24569,10.5)
df <- as.data.frame(mat)
colnames(df) <- c(" ","Black-Scholes","Heston","Market Quote")
print(df)

##                Black-Scholes   Heston Market Quote
## 1 In-the-Money      184.8947 189.0168        192.4
## 2 Out-of-Money             0 11.24569         10.5

The Black-Scholes results seem to be inaccurate. However, the options provided by APPLE Inc.(AAPL) are American type, the actual trading option price must be higher than or equal to the results computed by Black-Scholes equation(since APPL pays out dividends, the call options may be early exercised prior to the exdividend dates), where the formula can only deal with European option pricing.

Since the stochastic volatility models are calibrated to satisfy the implied volatility, they capture the effect- forward volatility. As a consequence, The Heston Model offers more precise prediction in option pricing. “Smile dynamics” is also captured more successfully in stochastic volatility models. Therefore, while pricing options, stochastic volatility models are preferred to Black-Scholes Model. However, when it comes to pricing the simple derivatives which only depend on the terminal distribution of the underlying, i.e. vanilla European options, Black-Scholes Model can provide a “good enough” result with a faster and easier way.

viii. Reference

Gatheral, J. and Taleb, N., 2011. The Volatility Surface. Hoboken: John Wiley & Sons, Inc.

Heston, S., 1993. A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. Review of Financial Studies, 6(2), pp.327-343.[Accessed 13 May 2020].

Yang, Yuan, 2013. “Valuing a European option with the Heston model.” Thesis. Rochester Institute of Technology.[Accessed 13 May 2020].