##' Generalized Black Scholes model for pricing vanilla European options
##'
##' Compute values of call and put options as well as the Greeks -
##' the sensitivities of the option price to various input arguments using the
##' Generalized Black Scholes model. "Generalized" means that the asset can
##' have a continuous dividend yield.
##'
##' The Generalized Black Scholes formula for call options is \cr
##' \eqn{e^{-r t} (s \; e^{g t} \; Nd1 - X \; Nd2)}{exp(-r * t) * (s * exp(g * t) * Nd1 - X * Nd2)} \cr
##' where \cr
##' \eqn{g = r - div\_yield} \cr
##' \eqn{Nd1 = N(d1)} and \eqn{Nd2 = N(d2)} \cr
##' \eqn{d1 = \frac{log(s / X) + (g + Sigma^2/ 2) t}{Sigma \sqrt{t}}}{d1 = (log(s / X) + (g + Sigma^2/ 2) * t) / (Sigma * sqrt(t))} \cr
##' \eqn{d2 = d1 - Sigma \sqrt{t}}{d2 = d1 - Sigma * sqrt(t)} \cr
##' N denotes the normal CDF (\code{\link{pnorm}})\cr
##' For put options, the formula is \cr
##' \eqn{e^{-r t} (-s \; e^{g t} \; Nminusd1 + X \; Nminusd2)}{exp(-r * t) * (-s * exp(g * t) * Nminusd1 + X * Nminusd2)}\cr
##' where \cr
##' \eqn{Nminusd1 = N(-d1)} and \eqn{Nminusd2 = N(-d2)} \cr
##'
##' @param s the spot price of the asset (the stock price for options on stocks)
##' @param X the exercise or strike price of the option
##' @param r the continuously compounded rate of interest in decimal (0.10 or 10e-2 for 10\%)
##' (use \code{\link{equiv.rate}} to convert to a continuously compounded rate)
##' @param Sigma the volatility of the asset price in decimal (0.20 or 20e-2 for 20\%)
##' @param t the maturity of the option in years
##' @param div_yield the continuously compounded dividend yield (0.05 or 5e-2 for 5\%)
##' (use \code{\link{equiv.rate}} to convert to a continuously compounded rate)
##' @return A list of the following elements
##' \item{call}{the value of a call option}
##' \item{put}{the value of a put option}
##' \item{Greeks}{a list of the following elements}
##' \item{Greeks$callDelta}{the delta of a call option - the sensitivity to the spot price of the asset}
##' \item{Greeks$putDelta}{the delta of a put option - the sensitivity to the spot price of the asset}
##' \item{Greeks$callTheta}{the theta of a call option - the time decay of the option value
##' with passage of time. Note that time is measured in years. To find a daily theta divided by 365.}
##' \item{Greeks$putTheta}{the theta of a put option}
##' \item{Greeks$Gamma}{the gamma of a call or put option - the second derivative with respect to the spot price
##' or the sensitivity of delta to the spot price}
##' \item{Greeks$Vega}{the vega of a call or put option - the sensitivity to the volatility}
##' \item{Greeks$callRho}{the rho of a call option - the sensitivity to the interest rate}
##' \item{Greeks$putRho}{the rho of a put option - the sensitivity to the interest rate}
##' \item{extra}{a list of the following elements}
##' \item{extra$d1}{the d1 of the Generalized Black Scholes formula}
##' \item{extra$d2}{the d2 of the Generalized Black Scholes formula}
##' \item{extra$Nd1}{is \code{\link{pnorm}}(d1)}
##' \item{extra$Nd2}{is \code{\link{pnorm}}(d2)}
##' \item{extra$Nminusd1}{is \code{\link{pnorm}}(-d1)}
##' \item{extra$Nminusd2}{is \code{\link{pnorm}}(-d2)}
##' \item{extra$callProb}{the (risk neutral) probability that the call will be exercised = \code{Nd2}}
##' \item{extra$putProb}{the (risk neutral) probability that the put will be exercised = \code{Nminusd2}}
##'
##' @export
##' @importFrom stats pnorm dnorm
GenBS <- function(s, X, r, Sigma, t, div_yield = 0){
g = r - div_yield ## the growth rate of the asset price
d1 <- (log(s / X) + (g + Sigma*Sigma / 2) * t) / (Sigma * sqrt(t))
## the division above can be 0/0=NaN if s=X and either t=0 or g=Sigma=0
## we set d1 to 0 in this case since numerator is of higher order in t and Sigma
d1 <- ifelse(is.nan(d1), 0, d1)
d2 <- d1 - Sigma * sqrt(t)
Nd1 <- pnorm(d1)
Nd2 <- pnorm(d2)
Nminusd1 <- pnorm(-d1)
Nminusd2 <- pnorm(-d2)
## Black Scholes call price
call <- exp(-r * t) * (s * exp(g * t) * Nd1 - X * Nd2)
put <- exp(-r * t) * (-s * exp(g * t) * Nminusd1 + X * Nminusd2)
callDelta <- exp(-div_yield*t)*Nd1
putDelta <- -exp(-div_yield*t)*Nminusd1
a <- -s * dnorm(d1) * Sigma * exp(-div_yield*t) / (2 * sqrt(t))
## if t is 0, the above is 0/0 = NaN but dnorm(d1) term goes to 0 faster
## so we set the result to 0 in this case
a <- ifelse(is.nan(a), 0, a)
b = r * X * exp(-r * t) * Nd2
c = div_yield * s * Nd1* exp(-div_yield*t)
callTheta <- a - b + c
b = r * X * exp(-r * t) * Nminusd2
c = div_yield * s * Nminusd1 * exp(-div_yield*t)
putTheta <- a + b - c
Gamma <- dnorm(d1) * exp(-div_yield*t)/ (s * Sigma * sqrt(t))
## if t is 0, the above is 0/0 = NaN but dnorm(d1) term goes to 0 faster
## so we set the result to 0 in this case
Gamma <- ifelse(is.nan(Gamma), 0, Gamma)
Vega <- s * sqrt(t) * dnorm(d1) * exp(-div_yield*t)
callRho <- X * t * exp(-r * t) * Nd2
putRho <- -X * t * exp(-r * t) * Nminusd2
callProb <- Nd2
putProb <- Nminusd2
Greeks <- list(callDelta=callDelta, putDelta=putDelta, callTheta=callTheta,
putTheta=putTheta, Gamma=Gamma, Vega=Vega, callRho=callRho, putRho=putRho)
extra <- list(d1=d1, d2=d2, Nd1=Nd1, Nd2=Nd2, Nminusd1=Nminusd1,
Nminusd2=Nminusd2, callProb=callProb, putProb=putProb)
list(call=call, put=put, Greeks=Greeks, extra=extra)
}
##' Generalized Black Scholes model implied volatility
##'
##' Find implied volatility given the option price using the generalized Black Scholes model.
##' "Generalized" means that the asset can have a continuous dividend yield.
##'
##' \code{GenBSImplied} calls \code{\link{newton.raphson.root}} and
##' if that fails \code{\link{uniroot}}
##'
##' @param price the price of the option
##' @param PutOpt \code{TRUE} for put options, \code{FALSE} for call options
##' @param toler passed on to \code{\link{newton.raphson.root}}
##' The implied volatility is regarded as correct if the solver is able to
##' match the option price to within less than \code{toler}. Otherwise the function returns \code{NA}
##' @param max.iter passed on to \code{\link{newton.raphson.root}}
##' @param convergence passed on to \code{\link{newton.raphson.root}}
##' @inheritParams GenBS
##'
##' @export
##' @importFrom stats uniroot
GenBSImplied <- function(s, X, r, price, t, div_yield, PutOpt=FALSE,
toler=1e-6, max.iter=100, convergence=1e-8){
## discount the exercise price to eliminate r
X = X * exp(-r * t)
## adjust the spot price to eliminate div_yield
s = s * exp(-div_yield * t)
## use put call parity to convert put option into call option
if (PutOpt) price = price + (s - X)
SminusX = s - X
SplusX = s + X
if (price < SminusX || price < 0 || price > s){
warning("Implied volatility is undefined because price is outside theoretical bounds")
return(NA)
}
if (price == SminusX || price == 0)
## if price equals intrinsic value, volatility is zero
return(0)
if (X == 0 && price != s) {
## if x is 0, option price must equal stock price
warning("Implied volatility is undefined because price is outside theoretical bounds")
return(NA)
}
## we use a Taylor series approximation for an initial guess
guess <- GenBSImpliedGuess(price, SminusX, SplusX, t)
## since price exceeds intrinsic value, 0 is a lower bound for the volatility
## we now seek an upper bound by doubling 100% volatility until the price is exceeded
upper <- 100e-2
while (GenBS(s, X, 0, upper, t, 0)$call < price && upper < .Machine$integer.max) upper <- upper * 2
if (upper > .Machine$integer.max) upper <- Inf
f <- function(sigma){
temp <- GenBS(s, X, 0, sigma, t, 0)
list(value = temp$call - price, gradient = temp$Greeks$Vega)
}
res <- newton.raphson.root(f=f, guess=guess, lower=0, upper=upper, max.iter=max.iter,
toler=toler, convergence=convergence)
if (! is.na(res))
return (res)
## if newton.raphson.root fails, we try bisection
res2 <- uniroot(function(x){f(x)$value}, c(0, upper), tol=convergence)
if (abs(res2$f.root > toler)){
warning("Error finding implied volatility")
return (NA) # uniroot also failed so return NA
}
warning("Found implied volatility using uniroot after newton.raphson.root failed")
## Above warning is to prevent confusion from warning message in newton.raphson.root
return (res2$root)
}
GenBSImpliedGuess <- function(price, SminusX, SplusX, t){
## We use a Taylor series approximation similar to
## Brenner, H. and Subrahmanyam, M. G. (1994), "A simple approximation to option valuation
## and hedging in the Black Scholes model", Financial Analysts Journal, Mar-Apr 1994, 25-28
## Assume a call option and that r has been eliminated by discounting the exercise price
## Approximate log(s/x) as 2(s-x)/(s+x) = SminusX(s-x)/H
## where SminusX = s-x and H=(s+x)/2.
## We then approximate the normal integral by a Taylor series
## if price is the call price, this gives the approximation
## price/H = SigmaRootT(1/Root2Pi +SminusX/(2*H*SigmaRootT)
## +SminusX^2/(2*H^2*SigmaRootT^2*Root2Pi)
## Letting z=SigmaRootT*H/Root2Pi, this yields a quadratic equation for z:
## z^2 - (price-SminusX/2)*z +SminusX^2/(4*Pi) = 0
## The solution is
## z = (price - SminusX/2)/2 + Radical
## where Radical is the square root of (price-SminusX/2)^2 - SminusX^2/Pi
## The linear approximation is obtained by dropping the constant term in the quadratic
## to give the linear equation:
## z = (price-SminusX/2)
root2pi = sqrt(2 * pi)
h = 0.5 * SplusX
temp = (price - 0.5 * SminusX)
radical = temp*temp - SminusX*SminusX / pi
if (radical < 0) # Try Linear Approximation
sigmaRootT = (root2pi / h) * temp
else{ # Try Quadratic Approximation
radical = sqrt(radical)
sigmaRootT = (root2pi / h) * (temp / 2 + radical)
}
sigmaRootT / sqrt(t)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.