Skip to content

Instantly share code, notes, and snippets.

@wilsonfreitas

wilsonfreitas/vol.R

Created March 2, 2018 22:35
Show Gist options
  • Save wilsonfreitas/b1532177263d1b60581227898aab6b63 to your computer and use it in GitHub Desktop.
Save wilsonfreitas/b1532177263d1b60581227898aab6b63 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
vol.R
#' Computes option price via Black-76 formula.
#'
#' Black-76 formula implementation, pricing european options.
#'
#' @param type 'call' for call option, any other value for put.
#' @param future current future price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param sigma volatility of the future price.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with price of european options on futures, according to Black-76 model.
#'
#' @seealso
#' \code{\link{bsmprice}} for pricing options on stocks, with
#' Black-Scholes-Merton model.
#'
#' @examples
#' blackprice(c('put', 'call'), c(50, 49, 48, 47), 45, 0.25, 0.13, 0.3)
#' blackprice('p', 100:500/10, 45, 0.25, 0.13, 0.17)
#' blackprice('call', 50, 100:900/10, 2, 0.1, 0.18)
#'
#' @export
blackprice <- function(type, future, strike, time, rate, sigma) {
df <- data.frame(type, future, strike, time, rate, sigma)
with(df, {
d1 <- (log(future / strike) + (sigma * sigma/2) * time) / (sigma * sqrt(time))
d2 <- d1 - sigma * sqrt(time)
ifelse(tolower(as.character(type)) == "call", {
future * exp((- rate) * time) * pnorm(d1) - strike * exp(- rate * time) * pnorm(d2)
}, {
strike * exp(- rate * time) * pnorm(- d2) - future * exp((rate) * time) * pnorm(-d1)
})
})
}
#' Computes the desired greek letter from Black-76 model.
#'
#' Greeks formula implementation, from Black-76 model.
#'
#' @param greeks desired greek letter: 'delta', 'gamma', 'vega', 'rho' or 'theta'.
#' @param type 'call' for call option, any other value for put.
#' @param future current underlying forward price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param sigma volatility of the stock price.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with desired greeks for input data.
#'
#' @examples
#' blackgreeks('delta', 'put', 20, 30, 0:2000/1000, 0.15, 0.25)
#'
#' @export
blackgreeks <- function(greeks=c('delta', 'gamma', 'vega', 'rho', 'theta'),
type, future, strike, time, rate, sigma) {
df <- data.frame(type, future, strike, time, rate, sigma)
greeks <- match.arg(greeks)
with(df, {
if(greeks=='delta') blackdelta(type, future, strike, time, rate, sigma)
if(greeks=='gamma') blackgamma(type, future, strike, time, rate, sigma)
if(greeks=='vega') blackvega(type, future, strike, time, rate, sigma)
if(greeks=='rho') blackrho(type, future, strike, time, rate, sigma)
if(greeks=='theta') blacktheta(type, future, strike, time, rate, sigma)
})
}
#' @rdname blackgreeks
#'
#' @section blackdelta:
#' Delta measures the rate of change of the theoretical option value with respect to changes in the underlying asset's price.
#'
#' @seealso
#' \code{\link{bsmdelta}} for delta from Black-Scholes-Merton model.
#'
#' @examples
#' blackdelta(c('call', 'put'), c(50, 49, 48, 47), c(45, 40), 0.25, 0.13, 0.2)
#' blackdelta('call', 50, 45, 0.25, 0.13, 1:1000/1000)
#' blackdelta('put', 20, 30, 0:2000/1000, 0.15, 0.25)
#'
#' @export
blackdelta <- function(type, future, strike, time, rate, sigma){
d1 <- (log(future / strike) + (sigma * sigma / 2) * time) / (sigma * sqrt(time))
if (tolower(as.character(type)) == "call"){
result <- exp((- rate) * time) * pnorm(d1)
} else {
result <- exp((- rate) * time) * (pnorm(d1) - 1)
}
result
}
#' @rdname blackgreeks
#'
#' @section blackvega:
#' Vega measures sensitivity of volatility.
#'
#' @seealso
#' \code{\link{bsmvega}} for vega from Black-Scholes-Merton model.
#'
#' @examples
#' blackvega(c('call', 'put'), c(50, 49, 48, 47), c(45, 40), 0.25, 0.13, 0.2)
#' blackvega('call', 50, 45, 0.25, 0.13, 1:1000/1000)
#' blackvega('put', 20, 30, 0:2000/1000, 0.15, 0.25)
#'
#' @export
blackvega <- function(type, future, strike, time, rate, sigma){
d1 <- (log(future / strike) + (sigma * sigma / 2) * time) / (sigma * sqrt(time))
result <- future * exp((- rate) * time) * dnorm(d1) * sqrt(time)
result
}
#' @rdname blackgreeks
#'
#' @section blackgamma:
#' Gamma measures the rate of change in the delta with respect to changes in the underlying price.
#'
#' @seealso
#' \code{\link{bsmgamma}} for gamma from Black-Scholes-Merton model.
#'
#' @examples
#' blackgamma(c('call', 'put'), c(50, 49, 48, 47), c(45, 40), 0.25, 0.13, 0.2)
#' blackgamma('call', 50, 45, 0.25, 0.13, 1:1000/1000)
#' blackgamma('put', 20, 30, 0:2000/1000, 0.15, 0.25)
#'
#' @export
blackgamma <- function(type, future, strike, time, rate, sigma){
d1 <- (log(future / strike) + (sigma * sigma / 2) * time) / (sigma * sqrt(time))
result <- exp((- rate) * time) * dnorm(d1)/(future * sigma * sqrt(time))
result
}
#' @rdname blackgreeks
#'
#' @section blacktheta:
#' Theta measures the sensitivity of the value of the derivative to the passage of time.
#'
#' @seealso
#' \code{\link{bsmtheta}} for theta from Black-Scholes-Merton model.
#'
#' @examples
#' blacktheta(c('call', 'put'), c(50, 49, 48, 47), c(45, 40), 0.25, 0.13, 0.2)
#' blacktheta('call', 50, 45, 0.25, 0.13, 1:1000/1000)
#' blacktheta('put', 20, 30, 0:2000/1000, 0.15, 0.25)
#'
#' @export
blacktheta <- function(type, future, strike, time, rate, sigma){
d1 <- (log(future / strike) + (sigma * sigma / 2) * time) / (sigma * sqrt(time))
d2 <- d1 - sigma * sqrt(time)
Theta1 <- -(future * exp((- rate) * time) * dnorm(d1) * sigma) / (2 * sqrt(time))
if (tolower(as.character(type)) == "call"){
result <- Theta1 - (- rate) * future * exp((- rate) * time) * pnorm(+ d1) - rate * strike * exp(- rate * time) * pnorm(+ d2)
} else {
result <- Theta1 + (- rate) * future * exp((- rate) * time) * pnorm(- d1) + rate * strike * exp(- rate * time) * pnorm(- d2)
}
result
}
#' @rdname blackgreeks
#'
#' @section blackrho:
#' Rho measures sensitivity to the interest rate: it is the derivative of the option value with respect to the risk free interest rate.
#'
#' @seealso
#' \code{\link{bsmrho}} for rho from Black-Scholes-Merton model.
#'
#' @examples
#' blackrho(c('call', 'put'), c(50, 49, 48, 47), c(45, 40), 0.25, 0.13, 0.2)
#' blackrho('call', 50, 45, 0.25, 0.13, 1:1000/1000)
#' blackrho('put', 20, 30, 0:2000/1000, 0.15, 0.25)
#'
#' @export
blackrho <- function(type, future, strike, time, rate, sigma){
d1 <- (log(future / strike) + (sigma * sigma / 2) * time) / (sigma * sqrt(time))
d2 <- d1 - sigma * sqrt(time)
CallPut = bsmprice(type, future, strike, time, rate, rate, sigma)
result <- - time * CallPut
result
}
#' Computes d1 from Black-76 model.
#'
#' d1 formula implementation, from Black-76 model.
#'
#' @param future current future price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param sigma volatility of the stock price.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with d1.
#'
#' @seealso
#' \code{\link{bsmd1}} for d1 from Black-Scholes-Merton model.
#'
#' @examples
#' blackd1(c(50, 49, 48, 47), 45, 0.25, 0.13, 0.2)
#' blackd1(50, 45, 0.25, 0.13, 1:1000/1000)
#' blackd1(20, 30, 0:2000/1000, 0.15, 0.25)
#'
#' @export
blackd1 <- function(future, strike, time, rate, sigma) {
(log(future / strike) + (sigma * sigma / 2) * time) / (sigma * sqrt(time))
}
#' Computes a measure of moneyness for Black-76 model.
#'
#' Moneyness formula implementation, for Black-76 model.
#'
#' @param future current future price.
#' @param strike the strike price.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with moneyness.
#'
#' @seealso
#' \code{\link{bsmmoneyness}} for one measure for moneyness in
#' Black-Scholes-Merton model.
#'
#' @examples
#' blackmoneyness(2:51, 1:50)
#' blackmoneyness(10, 50:150/10)
#' blackmoneyness(0:600/10, 30)
#'
#' @export
blackmoneyness <- function(future, strike) {
1 - future / strike
}
blackstrike <- function(future, delta, sigma, t) {
inv <- qnorm(delta)
ee <- sigma*sqrt(t)
strike <- future / ( exp((inv - 0.5*ee) * ee))
strike
}
#' Computes option price via Black-Scholes-Merton formula.
#'
#' Black-Scholes-Merton formula implementation, pricing european options.
#'
#' @param type 'call' for call option, any other value for put.
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#' @param sigma volatility of the stock price.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with price of european options on stocks, according to
#' Black-Scholes-Merton model.
#'
#' @seealso
#' \code{\link{blackprice}} for pricing options on futures, with Black-76 model.
#'
#' @examples
#' bsmprice(c('call', 'put'), c(50, 49, 48, 47), 45, 0.25, 0.13, 0.01, 0.2)
#' bsmprice('call', 50, 45, 0.25, 0.13, 0.01, 1:1000/1000)
#' bsmprice('put', 20, 30, 0:2000/1000, 0.15, 0, 0.25)
#'
#' @export
bsmprice <- function(type, spot, strike, time, rate, yield, sigma) {
df <- data.frame(type, spot, strike, time, rate, yield, sigma)
with(df, {
d1 <- (log(spot / strike) + (rate - yield + sigma * sigma/2) * time) / (sigma * sqrt(time))
d2 <- d1 - sigma * sqrt(time)
ifelse(tolower(as.character(type)) == "call", {
spot * exp((- yield) * time) * pnorm(d1) - strike * exp(- rate * time) * pnorm(d2)
}, {
strike * exp(- rate * time) * pnorm(- d2) - spot * exp((- yield) * time) * pnorm(-d1)
})
})
}
#' Computes d1 from Black-Scholes-Merton model.
#'
#' d1 formula implementation, from Black-Scholes-Merton model.
#'
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#' @param sigma volatility of the stock price.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with d1.
#'
#' @seealso
#' \code{\link{blackd1}} for d1 from Black-76 model.
#'
#' @examples
#' bsmd1(c(50, 49, 48, 47), 45, 0.25, 0.13, 0.01, 0.2)
#' bsmd1(50, 45, 0.25, 0.13, 0.01, 1:1000/1000)
#' bsmd1(20, 30, 0:2000/1000, 0.15, 0, 0.25)
#'
#' @export
bsmd1 <- function(spot, strike, time, rate, yield, sigma) {
(log(spot / strike) + (rate - yield + sigma * sigma / 2) * time) / (sigma * sqrt(time))
}
#' Computes the desired greek letter from Black-Scholes-Merton model.
#'
#' Greeks formula implementation, from Black-Scholes-Merton model.
#'
#' @param greeks desired greek letter: 'delta', 'gamma', 'vega', 'rho' or 'theta'.
#' @param type 'call' for call option, any other value for put.
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#' @param sigma volatility of the stock price.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with desired greeks for input data.
#'
#' @examples
#' bsmgreeks('delta', 'put', 20, 30, 0:2000/1000, 0.15, 0, 0.25)
#'
#' @export
bsmgreeks <- function(greeks=c('delta', 'gamma', 'vega', 'rho', 'theta'),
type, spot, strike, time, rate, yield, sigma) {
df <- data.frame(type, spot, strike, time, rate, yield, sigma)
greeks <- match.arg(greeks)
with(df, {
if(greeks=='delta') bsmdelta(type, spot, strike, time, rate, yield, sigma)
if(greeks=='gamma') bsmgamma(type, spot, strike, time, rate, yield, sigma)
if(greeks=='vega') bsmvega(type, spot, strike, time, rate, yield, sigma)
if(greeks=='rho') bsmrho(type, spot, strike, time, rate, yield, sigma)
if(greeks=='theta') bsmtheta(type, spot, strike, time, rate, yield, sigma)
})
}
#' @rdname bsmgreeks
#'
#' @section bsmdelta:
#' Delta greek letter measures the rate of change of the theoretical option value with respect to changes in the underlying asset's price.
#'
#' @seealso
#' \code{\link{blackdelta}} for delta from Black-76 model.
#'
#' @examples
#' bsmdelta(c('call', 'put'), c(50, 49, 48, 47), 45, 0.25, 0.13, 0.01, 0.2)
#' bsmdelta('call', 50, 45, 0.25, 0.13, 0.01, 1:1000/1000)
#' bsmdelta('put', 20, 30, 0:2000/1000, 0.15, 0, 0.25)
#'
#' @export
bsmdelta <- function(type, spot, strike, time, rate, yield, sigma){
d1 <- (log(spot / strike) + (rate - yield + sigma * sigma / 2) * time) / (sigma * sqrt(time))
if (tolower(as.character(type)) == "call"){
result <- exp((- yield) * time) * pnorm(d1)
} else {
result <- exp((- yield) * time) * (pnorm(d1) - 1)
}
result
}
#' @rdname bsmgreeks
#'
#' @section bsmvega:
#' Vega greek letter measures sensitivity of volatility.
#'
#' @seealso
#' \code{\link{blackvega}} for vega from Black-76 model.
#'
#' @examples
#' bsmvega(c('call', 'put'), c(50, 49, 48, 47), 45, 0.25, 0.13, 0.01, 0.2)
#' bsmvega('call', 50, 45, 0.25, 0.13, 0.01, 1:1000/1000)
#' bsmvega('put', 20, 30, 0:2000/1000, 0.15, 0, 0.25)
#'
#' @export
bsmvega <- function(type, spot, strike, time, rate, yield, sigma){
d1 <- (log(spot / strike) + (rate - yield + sigma * sigma / 2) * time) / (sigma * sqrt(time))
result <- spot * exp((- yield) * time) * dnorm(d1) * sqrt(time)
result
}
#' @rdname bsmgreeks
#'
#' @section bsmgamma:
#' Gamma greek letter measures the rate of change in the delta with respect to changes in the underlying price.
#'
#' @seealso
#' \code{\link{blackgamma}} for gamma from Black-76 model.
#'
#' @examples
#' bsmgamma(c('call', 'put'), c(50, 49, 48, 47), 45, 0.25, 0.13, 0.01, 0.2)
#' bsmgamma('call', 50, 45, 0.25, 0.13, 0.01, 1:1000/1000)
#' bsmgamma('put', 20, 30, 0:2000/1000, 0.15, 0, 0.25)
#'
#' @export
bsmgamma <- function(type, spot, strike, time, rate, yield, sigma){
d1 <- (log(spot / strike) + (rate - yield + sigma * sigma / 2) * time) / (sigma * sqrt(time))
result <- exp((- yield) * time) * dnorm(d1)/(spot * sigma * sqrt(time))
result
}
#' @rdname bsmgreeks
#'
#' @section bsmtheta:
#' Theta measures the sensitivity of the value of the derivative to the passage of time.
#
#' @seealso
#' \code{\link{blacktheta}} for theta from Black-76 model.
#'
#' @examples
#' bsmtheta(c('call', 'put'), c(50, 49, 48, 47), 45, 0.25, 0.13, 0.01, 0.2)
#' bsmtheta('call', 50, 45, 0.25, 0.13, 0.01, 1:1000/1000)
#' bsmtheta('put', 20, 30, 0:2000/1000, 0.15, 0, 0.25)
#'
#' @export
bsmtheta <- function(type, spot, strike, time, rate, yield, sigma){
d1 <- (log(spot / strike) + (rate - yield + sigma * sigma / 2) * time) / (sigma * sqrt(time))
d2 <- d1 - sigma * sqrt(time)
Theta1 <- -(spot * exp((- yield) * time) * dnorm(d1) * sigma) / (2 * sqrt(time))
if (tolower(as.character(type)) == "call"){
result <- Theta1 - (- yield) * spot * exp((- yield) * time) * pnorm(+ d1) -
rate * strike * exp(- rate * time) * pnorm(+ d2)
} else {
result <- Theta1 + (- yield) * spot * exp((- yield) * time) * pnorm(- d1) +
rate * strike * exp(- rate * time) * pnorm(- d2)
}
result
}
#' @rdname bsmgreeks
#'
#' @section bsmrho:
#' Rho measures sensitivity to the interest rate: it is the derivative of the option value with respect to the risk free interest rate.
#'
#' @seealso
#' \code{\link{blackrho}} for rho from Black-76 model.
#'
#' @examples
#' bsmrho(c('call', 'put'), c(50, 49, 48, 47), 45, 0.25, 0.13, 0.01, 0.2)
#' bsmrho('call', 50, 45, 0.25, 0.13, 0.01, 1:1000/1000)
#' bsmrho('put', 20, 30, 0:2000/1000, 0.15, 0, 0.25)
#'
#' @export
bsmrho <- function(type, spot, strike, time, rate, yield, sigma){
d1 <- (log(spot / strike) + (rate - yield + sigma * sigma / 2) * time) / (sigma * sqrt(time))
d2 <- d1 - sigma * sqrt(time)
CallPut = bsmprice(type, spot, strike, time, rate, yield, sigma)
if (tolower(as.character(type)) == "call") {
if (rate - yield != 0) {
result <- time * strike * exp(- rate * time) * pnorm(d2)
}
else {
result <- - time * CallPut
}
} else {
if (rate - yield != 0) {
result <- - time * strike * exp(- rate * time) * pnorm(- d2)
}
else {
result <- - time * CallPut
}
}
result
}
#' Computes a measure of moneyness for Black-Scholes-Merton model.
#'
#' Moneyness formula implementation, for Black-Scholes-Merton model.
#'
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with moneyness.
#'
#' @seealso
#' \code{\link{blackmoneyness}} for one measure for moneyness in Black-76 model.
#'
#' @examples
#' bsmmoneyness(0:49, 1:50, 3, 0.1, 0)
#' bsmmoneyness(10, 10, 1, 0.12, 0:1000/1000)
#' bsmmoneyness(10, 10, 1, 0:1000/1000, 0)
#'
#' @export
bsmmoneyness <- function(spot, strike, time, rate, yield) {
1 - spot * exp( - yield * time) / (strike * exp( - rate * time))
}
#' Computes DI future option price via Black-Scholes formula.
#'
#' Black-Scholes formula implementation, pricing DI future options.
#'
#' @param type 'call' for call option, 'put' for put option. May be abbreviated.
#' @param short_price unit price, based on the date of option expiration.
#' @param long_price unit price, based on the date of underlying asset expiration.
#' @param strike the strike price.
#' @param short_t_cur time to option expiration in years, based in current days.
#' @param long_t_cur time to underlying asset expiration in years, based in
#' current days.
#' @param short_t_biz time to option expiration in years, based in working days.
#' @param long_t_biz time to underlying asset expiration in years, based in
#' working days.
#' @param sigma volatility of the underlying asset price.
#'
#' @return
#' DI future option price.
#'
#' @seealso
#' \code{\link{bsmprice}} for pricing options on stocks, with
#' Black-Scholes-Merton model. \\
#' \code{\link{blackprice}} for pricing options
#' on futures, with Black-76 model.
#'
#' @examples
#' blackDI('call', 110000, 100000, 0.14, 1, 2, 1, 2, 0.15)
#' blackDI('p', 100000, 90000, 0.08, 1, 2, 1, 2, 0.15)
#'
#' @export
blackDI <- function(type = c("call", "put"), short_price, long_price, strike,
short_t_cur, long_t_cur, short_t_biz, long_t_biz, sigma) {
type <- match.arg(type)
df <- data.frame(short_price, long_price, strike,
short_t_cur, long_t_cur, short_t_biz, long_t_biz, sigma)
df <- within(df, {
Kast <- ((1 + strike) ^ (long_t_biz - short_t_biz) - 1) * 1 / (long_t_cur - short_t_cur)
Szero <- (short_price / long_price - 1) * 1 / (long_t_cur - short_t_cur)
delta <- (long_price * (long_t_cur - short_t_cur)) / (1 + Kast * (long_t_cur - short_t_cur))
d1 <- (log(Szero / Kast) + (sigma * sigma) / 2 * short_t_cur) / (sigma * (short_t_cur) ^ 0.5)
d2 <- d1 - sigma * (short_t_cur) ^ 0.5
})
with(df,
ifelse(type == "call",
delta * (Szero * pnorm(d1) - Kast * pnorm(d2)),
delta * ( - Szero * pnorm( - d1) + Kast * pnorm( - d2))
))
}
#' @export
DELTAS <- c(1, 5, 10, 25, 35, 37, 50, 65, 67, 75, 90, 95, 99)
# Default SABR initial parameters
#' @export
ALPHA <- 0.001
#' @export
NU <- 0.001
#' @export
RHO <- 0.3
#' @export
BETA <- 0
# Default SVI initial parameters
#' @export
SVI.A <- 1e-12
#' @export
SVI.B <- 1e-8
#' @export
SVI.SIGMA <- 1e-8
#' @export
SVI.RHO <- -0.8
#' @export
SVI.M <- 1e-8
#' @export
SVI.ERR <- 0.02
# Default corrado-su initial parameters
#' @export
CS.SKEW <- 1e-8
#' @export
CS.KURT <- 1
#' @export
CS.SIGMA <- 1e-8
#' Computes w from modified Corrado-Su model.
#'
#' w formula implementation from modified Corrado-Su model.
#'
#' @param sigma volatility of the stock price.
#' @param time time to option expiration in years.
#' @param mu3 skewness.
#' @param mu4 kurtosis.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with w.
#'
#' @seealso
#' \code{\link{csmd}} for d from modified Corrado-Su model.\\
#' \code{\link{csmq3}} for Q_3 from modified Corrado-Su model.\\
#' \code{\link{csmq4}} for Q_4 from modified Corrado-Su model.\\
#' \code{\link{csmprice}} for modified Corrado-Su option pricing formula.
#'
#' @examples
#' csmw(0.25, 1.5, 0, -100:100/100)
#' csmw(0.25, 1.5, -100:100/100, 3)
#' csmw(0:100/100, 1.5, 0, 3)
#'
#' @export
csmw <- function(sigma, time, mu3, mu4) {
(mu3 / 6) * (sigma ^ 3) * (sqrt(time ^ 3)) + (mu4 / 24) * (sigma ^ 4) * (time ^ 2)
}
#' Computes d from modified Corrado-Su model.
#'
#' d formula implementation from modified Corrado-Su model.
#'
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#' @param sigma volatility of the stock price.
#' @param w constant from modified Corrado-Su model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with d.
#'
#' @seealso
#' \code{\link{csmw}} for w from modified Corrado-Su model.\\
#' \code{\link{csmq3}} for Q_3 from modified Corrado-Su model.\\
#' \code{\link{csmq4}} for Q_4 from modified Corrado-Su model.\\
#' \code{\link{csmprice}} for modified Corrado-Su option pricing formula.
#'
#' @examples
#' csmd(14, 120:160/10, 0.5, 0.15, 0.1, 0.18, csmw(0.18, 0.5, 0, 3))
#' csmd(14, 13, 0.5, 0.15, 0.1, 0.18, csmw(0.18, 0.5, -1000:1000/100, 3))
#' csmd(14, 13, 0.5, 0.15, 0.1, 0.18, csmw(0.18, 0.5, 0, -1000:1000/100))
#'
#' @export
csmd <- function(spot, strike, time, rate, yield, sigma, w) {
(log(spot / strike) + (rate - yield + 0.5 * sigma ^ 2) * time - log(1 + w)) / (sigma * sqrt(time))
}
#' Computes Q_3 from modified Corrado-Su model.
#'
#' Q_3 formula implementation from modified Corrado-Su model.
#'
#' @param spot current stock price.
#' @param sigma volatility of the stock price.
#' @param time time to option expiration in years.
#' @param d1 d value from modified Corrado-Su model.
#' @param w constant from modified Corrado-Su model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with Q_3.
#'
#' @seealso
#' \code{\link{csmw}} for w from modified Corrado-Su model.\\
#' \code{\link{csmd}} for d from modified Corrado-Su model.\\
#' \code{\link{csmq4}} for Q_4 from modified Corrado-Su model.\\
#' \code{\link{csmprice}} for modified Corrado-Su option pricing formula.
#'
#' @examples
#' spot <- 25
#' sigma <- 0.2
#' time <- 0.5
#' strike <- 26
#' rate <- 0.12
#' yield <- 0
#' skewness <- 0
#' kurtosis <- 3
#' w <- csmw(sigma, time, skewness, kurtosis)
#' d <- csmd(spot, strike, time, rate, yield, sigma, w)
#' csmq3(spot, sigma, time, d, w)
#' csmq3(200:320/10, sigma, time, d, w)
#'
#' @export
csmq3 <- function(spot, sigma, time, d1, w){
q3 <- spot * sigma * sqrt(time) / (6 * (1 + w))
q3 * (2 * sigma * sqrt(time) - d1) * dnorm(d1)
}
#' Computes Q_4 from modified Corrado-Su model.
#'
#' Q_4 formula implementation from modified Corrado-Su model.
#'
#' @param spot current stock price.
#' @param sigma volatility of the stock price.
#' @param time time to option expiration in years.
#' @param d1 d value from modified Corrado-Su model.
#' @param w constant from modified Corrado-Su model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with Q_4.
#'
#' @seealso
#' \code{\link{csmw}} for w from modified Corrado-Su model.\\
#' \code{\link{csmd}} for d from modified Corrado-Su model.\\
#' \code{\link{csmq3}} for Q_3 from modified Corrado-Su model.\\
#' \code{\link{csmprice}} for modified Corrado-Su option pricing formula.
#'
#' @examples
#' spot <- 25
#' sigma <- 0.2
#' time <- 0.5
#' strike <- 26
#' rate <- 0.12
#' yield <- 0
#' skewness <- 0
#' kurtosis <- 3
#' w <- csmw(sigma, time, skewness, kurtosis)
#' d <- csmd(spot, strike, time, rate, yield, sigma, w)
#' csmq4(spot, sigma, time, d, w)
#' csmq4(200:320/10, sigma, time, d, w)
#'
#' @export
csmq4 <- function(spot, sigma, time, d1, w){
q4 <- (d1 ^ 2 - 3 * d1 * sigma * sqrt(time) + 3 * sigma ^ 2 * time - 1) * dnorm(d1)
q4 <- q4 * spot * sigma * sqrt(time) / (24 * (1 + w))
return(q4)
}
#' Computes option price with modified Corrado-Su model.
#'
#' Implementation of modified Corrado-Su formula for option pricing.
#'
#' @param type 'call' for call option, any other value for put.
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#' @param sigma volatility of the stock price.
#' @param mu3 skewness.
#' @param mu4 kurtosis.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with option prices.
#'
#' @seealso
#' \code{\link{csmw}} for w from modified Corrado-Su model.\\
#' \code{\link{csmd}} for d from modified Corrado-Su model.\\
#' \code{\link{csmq3}} for Q_3 from modified Corrado-Su model.\\
#' \code{\link{csmq4}} for Q_4 from modified Corrado-Su model.
#'
#' @examples
#' spot <- 25
#' sigma <- 0.2
#' time <- 0.5
#' strike <- 26
#' rate <- 0.12
#' yield <- 0
#' skewness <- 1
#' kurtosis <- 4
#' csmprice(c('call', 'put'), spot, strike, time, rate, yield, sigma, 0, 3) ==
#' bsmprice(c('call', 'put'), spot, strike, time, rate, yield, sigma)
#' csmprice('put', spot, strike, time, rate, 0:100/100, sigma, skewness, kurtosis)
#'
#' @export
csmprice <- function(type, spot, strike, time, rate, yield, sigma, mu3, mu4) {
c_bs <- bsmprice('call', spot, strike, time, rate, yield, sigma)
w <- csmw(sigma, time, mu3, mu4)
d1 <- csmd(spot, strike, time, rate, yield, sigma, w)
q3 <- csmq3(spot, sigma, time, d1, w)
q4 <- csmq4(spot, sigma, time, d1, w)
callp <- c_bs + mu3 * q3 + (mu4 - 3) * q4
df <- data.frame(callp, type, spot, yield, time, rate, strike)
with(df, ifelse(tolower(as.character(type)) == 'call', callp, callp - spot * exp(- yield * time) + strike * exp(- rate * time)))
}
#' Computes Vega from modified Corrado-Su model.
#'
#' This function computes vega greek for an option based on modified Corrado-Su
#' pricing model.
#'
#' @param type 'call' for call option, any other value for put.
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#' @param sigma volatility of the stock price.
#' @param mu3 skewness.
#' @param mu4 kurtosis.
#'
#' @return
#' Vector with options Vega.
#'
#' @examples
#' csmvega('call', 14, 13, 0.1, 0.15, 0, 0.5, 0.18, 3.53)
#' csmvega('put', 14, 13, 0.1, 0.15, 0, 0.5, 0.18, 2.53)
#'
#' @export
csmvega <- function(type, spot, strike, time, rate, yield, sigma, mu3, mu4){
dsig <- 0.001
csmu <- csmprice(type, spot, strike, time, rate, yield, sigma + dsig, mu3, mu4)
csmd <- csmprice(type, spot, strike, time, rate, yield, sigma - dsig, mu3, mu4)
return((csmu - csmd) / (2 * dsig))
}
#' Computes Objective Function for Optimizing Modified Corrado-Su Parameters.
#'
#' This is an auxiliary function that provides the objective function used in
#' the optimization to obtain the Modified Corrado-Su parameters.
#'
#' @param par a 3-dimension vector with the parameters
#' @param spot current stock price.
#' @param strike the strike price.
#' @param rate the risk-free interest rate.
#' @param time time to option expiration in years.
#' @param y.data known prices for options.
#' @param sy.data known standard deviation for prices for options.
#'
#' @return
#' Normalized error between prices obtained from data and the ones obtained
#' from model.
#'
#' @export
f.optim.csm <- function(par, spot, strike, rate, time, y.data, sy.data){
sigma <- par[1]
mu3 <- par[2]
mu4 <- par[3]
yf <- csmprice(spot, strike, sigma, rate, time, mu3, mu4)
return(sum(((yf - y.data) / sy.data) ^ 2))
}
#' Computes Inequality Constraint for Optimizing Modified Corrado-Su Parameters.
#'
#' This is an auxiliary function that provides the inequality constraint used in
#' the optimization to obtain the Modified Corrado-Su parameters.
#'
#' @param par a 3-dimension vector with the parameters
#' @param spot current stock price.
#' @param strike the strike price.
#' @param rate the risk-free interest rate.
#' @param time time to option expiration in years.
#' @param y.data known prices for options.
#' @param sy.data known standard deviation for prices for options.
#'
#' @return
#' Constraint to be considered in the optimization.
#'
#' @export
g_ineq <- function(par, spot, strike, rate, time, y.data, sy.data){
mu3 <- par[2]
mu4 <- par[3]
return(abs(mu3) - skew_gram(mu4))
#tst <- abs(mu3)-skew_gram(mu4)
#if(tst > 0){
# print(paste(mu3, mu4, tst, sep=' '))
#}
#return(tst)
}
#' Computes Gradient of Inequality Constraint for Optimizing Modified Corrado-Su Parameters.
#'
#' This is an auxiliary function that provides the gradient for inequality constraint
#' used in the optimization to obtain the Modified Corrado-Su parameters.
#'
#' @param par a 3-dimension vector with the parameters
#' @param spot current stock price.
#' @param strike the strike price.
#' @param rate the risk-free interest rate.
#' @param time time to option expiration in years.
#' @param y.data known prices for options.
#' @param sy.data known standard deviation for prices for options.
#'
#' @return
#' Three-dimension vector with the gradient of the inequality constraint.
#'
#' @export
grad_ineq <- function(par, spot, strike, rate, time, y.data, sy.data){
return(c(0, 1, 1))
}
#' Computes Gradient of the Objective Function for Optimizing Modified Corrado-Su Parameters.
#'
#' This is an auxiliary function that provides the gradient of the objective
#' function used in the optimization to obtain the Modified Corrado-Su parameters.
#'
#' @param par a 3-dimension vector with the parameters
#' @param type 'call' for call option, any other value for put.
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#' @param y.data known prices for options.
#' @param sy.data known standard deviation for prices for options.
#'
#' @return
#' Matrix with the gradient components of the objective function.
#'
#' @export
grad <- function(par, type, spot, strike, time, rate, yield, y.data, sy.data){
#numeric vega for Corrado-Su mod
dcs.dsig <- csmvega(type, spot, strike, time, rate, yield,
par[1], par[2], par[3])
#other derivatives for the gradient
csmw <- csmw(par[1], time, par[2], par[3])
dmod <- csmd(spot, strike, 0, par[1], time, csmw)
q3 <- csmq3(spot, par[1], time, dmod, csmw)
q4 <- csmq4(spot, par[1], time, dmod, csmw)
premium <- csmprice(spot, strike, par[1], rate, time, par[2], par[3])
grad1 <- sum(2 * dcs.dsig * (premium - y.data) / (sy.data ^ 2))
grad2 <- sum(2 * q3 * (premium - y.data) / (sy.data ^ 2))
grad3 <- sum(2 * q4 * (premium - y.data) / (sy.data ^ 2))
return(cbind(grad1, grad2, grad3))
}
#' Computes Objective Function for Optimizing Bi-Modal Corrado-Su Parameters.
#'
#' This is an auxiliary function that provides the objective function used in
#' the optimization to obtain the Bi-Modal Corrado-Su parameters.
#'
#' @param par a 7-dimension vector with the parameters
#' @param spot current stock price.
#' @param strike the strike price.
#' @param rate the risk-free interest rate.
#' @param time time to option expiration in years.
#' @param y.data known prices for options.
#' @param sy.data known standard deviation for prices for options.
#'
#' @return
#' Normalized error between prices obtained from data and the ones obtained
#' from model.
#'
#' @export
f.optim.2csm <- function(par, spot, strike, rate, time, y.data, sy.data){
sigma1 <- par[1]
mu31 <- par[2]
mu41 <- par[3]
sigma2 <- par[4]
mu32 <- par[5]
mu42 <- par[6]
w <- par[7]
yf1 <- csmprice(spot, strike, sigma1, rate, time, mu31, mu41)
yf2 <- csmprice(spot, strike, sigma2, rate, time, mu32, mu42)
yf <- yf1 * w + yf2 * (1 - w)
return(sum(((yf - y.data) / sy.data) ^ 2))
}
#' @rdname f.optim.2csm
#' @export
csm.2opt <- function(par, spot, strike, rate, time){
sigma1 <- par[1]
mu31 <- par[2]
mu41 <- par[3]
sigma2 <- par[4]
mu32 <- par[5]
mu42 <- par[6]
w <- par[7]
yf1 <- csmprice(spot, strike, sigma1, rate, time, mu31, mu41)
yf2 <- csmprice(spot, strike, sigma2, rate, time, mu32, mu42)
return(yf1 * w + yf2 * (1 - w))
}
#' Computes Inequality Constraint for Optimizing Bi-Modal Corrado-Su Parameters.
#'
#' This is an auxiliary function that provides the inequality constraint used in
#' the optimization to obtain the Bi-Modal Corrado-Su parameters.
#'
#' @param par a 7-dimension vector with the parameters
#' @param spot current stock price.
#' @param strike the strike price.
#' @param rate the risk-free interest rate.
#' @param time time to option expiration in years.
#' @param y.data known prices for options.
#' @param sy.data known standard deviation for prices for options.
#'
#' @return
#' Constraint to be considered in the optimization.
#'
#' @export
g_ineq2 <- function(par, spot, strike, rate, time, y.data, sy.data){
mu31 <- par[2]
mu41 <- par[3]
mu32 <- par[5]
mu42 <- par[6]
tst1 <- abs(mu31) - skew_gram(mu41)
tst2 <- abs(mu32) - skew_gram(mu42)
if(tst1 > 0 || tst2 > 0){
#print('tst > 0')
return(1)
}else{
#print('tst < 0')
return(-1)
}
}
#' Computes Gradient of Inequality Constraint for Optimizing Bi-Modal Corrado-Su Parameters.
#'
#' This is an auxiliary function that provides the gradient for inequality constraint
#' used in the optimization to obtain the Bi-Modal Corrado-Su parameters.
#'
#' @param par a 7-dimension vector with the parameters
#' @param spot current stock price.
#' @param strike the strike price.
#' @param rate the risk-free interest rate.
#' @param time time to option expiration in years.
#' @param y.data known prices for options.
#' @param sy.data known standard deviation for prices for options.
#'
#' @return
#' Three-dimension vector with the gradient of the inequality constraint.
#'
#' @export
grad_ineq2 <- function(par, spot, strike, rate, time, y.data, sy.data){
return(c(0, 1, 1, 0, 1, 1, 0))
}
#' Computes Gradient of the Objective Function for Optimizing Bi-Modal Corrado-Su Parameters.
#'
#' This is an auxiliary function that provides the gradient of the objective
#' function used in the optimization to obtain the Bi-Modal Corrado-Su parameters.
#'
#' @param par a 7-dimension vector with the parameters
#' @param type 'call' for call option, any other value for put.
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#' @param y.data known prices for options.
#' @param sy.data known standard deviation for prices for options.
#'
#' @return
#' Matrix with the gradient components of the objective function.
#'
#' @export
grad2 <- function(par, type, spot, strike, time, rate, yield, y.data, sy.data){
#numeric vega for Corrado-Su mod
dcs.dsig1 <- csmvega(type, spot, strike, time, rate,
yield, par[1], par[2], par[3])
dcs.dsig2 <- csmvega(type, spot, strike, time, rate,
yield, par[4], par[5], par[6])
#other derivatives for the gradient
csmw1 <- csmw(par[1], time, par[2], par[3])
dmod1 <- csmd(spot, strike, 0, par[1], time, csmw1)
q31 <- csmq3(spot, par[1], time, dmod1, csmw1)
q41 <- csmq4(spot, par[1], time, dmod1, csmw1)
csmw2 <- csmw(par[4], time, par[5], par[6])
dmod2 <- csmd(spot, strike, 0, par[4], time, csmw2)
q32 <- csmq3(spot, par[4], time, dmod2, csmw2)
q42 <- csmq4(spot, par[4], time, dmod2, csmw2)
premium1 <- csmprice(spot, strike, par[1], rate, time, par[2], par[3])
premium2 <- csmprice(spot, strike, par[4], rate, time, par[5], par[6])
grad1 <- sum(2 * dcs.dsig1 * par[7] * (premium1 - y.data) / (sy.data ^ 2))
grad2 <- sum(2 * q31 * par[7] * (premium1 - y.data) / (sy.data ^ 2))
grad3 <- sum(2 * q41 * par[7] * (premium1 - y.data) / (sy.data ^ 2))
grad4 <- sum(2 * dcs.dsig2 * par[7] * (premium2 - y.data) / (sy.data ^ 2))
grad5 <- sum(2 * q32 * par[7] * (premium2 - y.data) / (sy.data ^ 2))
grad6 <- sum(2 * q42 * par[7] * (premium2 - y.data) / (sy.data ^ 2))
grad7 <- sum((premium1 - premium2) * (premium2 - y.data) / (sy.data ^ 2))
#print('grd')
#print(cbind(grad1, grad2, grad3, grad4, grad5, grad6, grad7))
return(cbind(grad1, grad2, grad3, grad4, grad5, grad6, grad7))
}
#' Modified Corrado-Su std. Deltas
#'
#' Finds and returns a dataframe with the strikes and implied
#' volatilites for the provided set of deltas or for the
#' default deltas.
#'
#' @param type Option type, "call" for calls or anything else for puts
#' @param spot underlying spot price at a given date
#' @param time time in years from a given date to a maturity date
#' @param rate interest free rate
#' @param yield dividend yield
#' @param sigma volatility of the stock price
#' @param mu3 skewness
#' @param mu4 kurtosis
#' @param deltas deltas for which strikes and vols should be calculated, 1-99. Defaults to usual
#' delta values.
#'
#' @return
#' Dataframe with columns delta, strike and impvol
#' @export
csm_deltas <- function(type, spot, time, rate, yield, sigma, mu3, mu4, deltas = DELTAS) {
ft <- spot * exp((rate - yield) * time)
eps <- 1e-8
strk_cs <- sapply(deltas/100, function(delta) {
uniroot(f=csm_strikes, interval=c(1, spot*10), tol=eps,
delta=delta, type=type, spot=spot, time=time,
rate=rate, yield=yield, sigma=sigma, mu3=mu3, mu4=mu4, ft=ft)$root
})
impvol_cs <-blackimpvol(csmprice(type, spot, strk_cs, time, rate, yield, sigma, mu3, mu4),
type, ft, strk_cs, time, rate)
data.frame(delta=deltas, strike=strk_cs, impvol=impvol_cs)
}
csm_strikes <- function(x, delta, type, spot, strike, time, rate, yield, sigma, mu3, mu4, ft) {
csp <- csmprice(type, spot, x, time, rate, yield, sigma, mu3, mu4)
if (csp < 0) return(999)
varx <- blackimpvol(csp, type, ft, x, time, rate)
x - blackstrike(ft, delta, varx, time)
}
#' Check Gram Charlier limits from parameters
#'
#' This function verifies if kurtosis and skewness of a distribution satisfy
#' Gram-Charlier Kurtosis-Skewness limits.
#'
#' @param k kurtosis.
#' @param s skewness.
#' @param verbose TRUE or FALSE for whether to show or not warnings in console.
#'
#' @return
#' TRUE or FALSE, in or out of the limits, respectively.
#'
#' @seealso
#' \code{\link{check_gram_charlier_limits_4_series}}, \code{\link{adjmoments}}.
#'
#' @examples
#' check_gram_charlier_limits(4, .5)
#' check_gram_charlier_limits(6, -.5)
#' check_gram_charlier_limits(3, 0)
#'
#' @export
check_gram_charlier_limits <- function(k, s, verbose=FALSE) {
data(gram_charlier_limits)
rk <- range(gram_charlier_limits$points[,1])
if ((k < rk[1] || k > rk[2]) && verbose)
warning('kurtosis exceeded: ', k)
rs <- max(abs(gram_charlier_limits$points[,2]))
if (s > rs && verbose)
warning('skewness exceeded: ', s)
x <- gram_charlier_limits$fun(k) > s
if (is.na(x)) FALSE else x
}
#' Check Gram Charlier limits from series
#'
#' This function verifies if kurtosis and skewness of a distribution satisfy
#' Gram-Charlier Kurtosis-Skewness limits.
#'
#' @param series return series to verify limits.
#' @param verbose TRUE or FALSE for whether to show or not warnings in console.
#'
#' @return
#' TRUE or FALSE, in or out of the limits, respectively.
#'
#' @seealso
#' \code{\link{check_gram_charlier_limits}}, \code{\link{adjmoments}}.
#'
#' @examples
#' check_gram_charlier_limits_4_series(rnorm(100000))
#'
#' @export
check_gram_charlier_limits_4_series <- function(series, verbose=FALSE) {
k <- kurtosis(series, method='moment')
s <- abs(skewness(series, method='moment'))
check_gram_charlier_limits(k, s, verbose)
}
#' Out of the Boundaries Skewness and Kurtosis Adjustment
#'
#' In case skewness and kurtosis fall outside Gram-Charlier limits, use this
#' function to adjust them so they fall inside the borders of the region.
#'
#' @param s skewness
#' @param k kurtosis
#'
#' @return
#' Dataframe with adjusted skewness and kurtosis.
#'
#' @seealso
#' \code{\link{check_gram_charlier_limits}}, \code{\link{adjmoments}}.
#'
#' @examples
#' adjmoments(-200:200/100, 8)
#' adjmoments(0.5, 10:90/10)
#'
#' @export
adjmoments <- function(s, k){
data <- data.frame(s, k)
if(nrow(data)==0) return()
data(gram_charlier_limits)
edges <- as.data.frame(gram_charlier_limits$points)
data <- within(data,{
adj.skewness <- NA
adj.kurtosis <- NA
})
na <- with(data, (is.na(s) | is.na(k)))
kurtok <- with(data, !na & k >= 3 & k <= 7)
data.kurtok <- data[kurtok,]
data.kurtok <- within(data.kurtok, {
adj.kurtosis <- k
skew.max <- approx(edges[,c(1,3)], xout = k)[[2]]
adj.skewness <- ifelse(abs(s) <= skew.max, s,
sign(s) * skew.max)
})
data[kurtok,] <- data.kurtok[,names(data)]
sub.edges <- list(subset(edges[order(edges[,3]),], k <= 5),
subset(edges[order(edges[,3]),], k > 5))
skewok <- with(data, !na & !kurtok & s >= -1 & s <= 1)
data.skewok <- data[skewok,]
data.skewok <- within(data.skewok, {
adj.skewness <- s
adj.kurtosis <- ifelse(k < 3,
approx(sub.edges[[1]][,c(3,1)], xout=abs(s), ties = "ordered")[[2]],
approx(sub.edges[[2]][,c(3,1)], xout=abs(s), ties = "ordered")[[2]])
})
data[skewok,] <- data.skewok[,names(data)]
nok <- with(data, !na & !kurtok & !skewok)
data.nok <- data[nok,]
data.nok <- within(data.nok, {
adj.kurtosis <- ifelse(k > 7, 7, 3)
adj.skewness <- sign(s) * 0.5
})
data[nok,] <- data.nok[,names(data)]
return(data[,c('adj.skewness','adj.kurtosis')])
}
#' Computes Q_3 from Corrado-Su model.
#'
#' Implementation of Q_3 formula from Corrado-Su model. Q_3 measures the effects
#' of nonnormal skewness on the option price given by the model.
#'
#' @param spot current stock price.
#' @param sigma volatility of the stock price.
#' @param time time to option expiration in years.
#' @param d1 auxiliary value, same from Black-Scholes-Merton model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with Q_3.
#'
#' @seealso
#' \code{\link{csq4}} for Q_4 formula.\\
#' \code{\link{csprice}} for Corrado-Su option pricing formula.
#'
#' @examples
#' csq3(10, .3, 1.5, 0:100/100)
#' spot <- 25
#' sigma <- 0.2
#' time <- 0.5
#' strike <- 26
#' rate <- 0.12
#' yield <- 0
#' d1 <- bsmd1(spot, strike, time, rate, yield, sigma)
#' csq3(spot, sigma, time, d1)
#'
#' @export
csq3 <- function(spot, sigma, time, d1) {
q3 <- spot * sigma * sqrt(time) / 6
q3 * ((2 * sigma * sqrt(time) - d1) * dnorm(d1) + sigma * sigma * time * pnorm(d1))
}
#' Computes Q_4 from Corrado-Su model.
#'
#' Implementation of Q_4 formula from Corrado-Su model. Q_4 measures the effects
#' of nonnormal kurtosis on the option price given by the model.
#'
#' @param spot current stock price.
#' @param sigma volatility of the stock price.
#' @param time time to option expiration in years.
#' @param d1 auxiliary value, same from Black-Scholes-Merton model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with Q_4.
#'
#' @seealso
#' \code{\link{csq3}} for Q_3 formula.\\
#' \code{\link{csprice}} for Corrado-Su option pricing formula.
#'
#' @examples
#' csq4(10, .3, 1.5, 0:100/100)
#' spot <- 25
#' sigma <- 0.2
#' time <- 0.5
#' strike <- 26
#' rate <- 0.12
#' yield <- 0
#' d1 <- bsmd1(spot, strike, time, rate, yield, sigma)
#' csq4(spot, sigma, time, d1)
#'
#' @export
csq4 <- function(spot, sigma, time, d1) {
q4 <- (d1 ^ 2 - 1 - 3 * sigma * sqrt(time) * (d1 - sigma * sqrt(time))) * dnorm(d1)
q4 <- q4 + (sigma * sigma * sigma * time ^ (1.5) * pnorm(d1))
q4 * spot * sigma * sqrt(time) / 24
}
#' Computes option price with Corrado-Su model.
#'
#' Implementation of Corrado-Su formula for option pricing.
#'
#' @param type 'call' for call option, any other value for put.
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#' @param sigma volatility of the stock price.
#' @param mu3 skewness.
#' @param mu4 kurtosis.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with option prices.
#'
#' @seealso
#' \code{\link{csq3}} for Q_3 formula from Corrado-Su model.\\
#' \code{\link{csq4}} for Q_4 formula from Corrado-Su model.
#'
#' @examples
#' spot <- 25
#' sigma <- 0.2
#' time <- 0.5
#' strike <- 26
#' rate <- 0.12
#' yield <- 0
#' csprice('call', spot, strike, time, rate, yield, sigma, 0, 3) ==
#' bsmprice('call', spot, strike, time, rate, yield, sigma)
#' csprice('call', spot, strike, time, rate, yield, sigma, 0:100/100, 3)
#' csprice('call', spot, strike, time, rate, yield, sigma, 0.5, 0:1000/100)
#'
#' @export
csprice <- function(type, spot, strike, time, rate, yield, sigma, mu3, mu4) {
c_bs <- bsmprice('call', spot, strike, time, rate, yield, sigma)
d1 <- bsmd1(spot, strike, time, rate, yield, sigma)
q3 <- csq3(spot, sigma, time, d1)
q4 <- csq4(spot, sigma, time, d1)
callp <- c_bs + mu3 * q3 + (mu4 - 3) * q4
df <- data.frame(callp, type, spot, yield, time, rate, strike)
with(df, ifelse(tolower(as.character(type)) == 'call', callp, callp - spot * exp(- yield * time) + strike * exp(- rate * time)))
}
#' Computes Q_3 from Corrado-Su model options on futures.
#'
#' Implementation of Q_3 formula from Corrado-Su model for futures. Q_3 measures
#' the effects of nonnormal skewness on the option price given by the model.
#'
#' @param future current future price.
#' @param sigma volatility of the future price.
#' @param time time to option expiration in years.
#' @param d1 auxiliary value, same from Black-76 model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with Q_3.
#'
#' @seealso
#' \code{\link{blackcsq4}} for Q_4 formula for futures.\\
#' \code{\link{blackcsprice}} for Corrado-Su option pricing formula for futures.
#'
#' @examples
#' blackcsq3(10, .3, 1.5, 0:100/100)
#' future <- 25
#' sigma <- 0.2
#' time <- 0.5
#' strike <- 26
#' rate <- 0.12
#' d1 <- blackd1(future, strike, time, rate, sigma)
#' blackcsq3(future, sigma, time, d1)
#'
#' @export
blackcsq3 <- function(future, sigma, time, d1) {
q3 <- future * sigma * sqrt(time) / 6
q3 * ((2 * sigma * sqrt(time) - d1) * dnorm(d1) + sigma * sigma * time * pnorm(d1))
}
#' Computes Q_4 from Corrado-Su model options on futures.
#'
#' Implementation of Q_4 formula from Corrado-Su model for futures. Q_4 measures
#' the effects of nonnormal skewness on the option price given by the model.
#'
#' @param future current future price.
#' @param sigma volatility of the future price.
#' @param time time to option expiration in years.
#' @param d1 auxiliary value, same from Black-76 model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with Q_4.
#'
#' @seealso
#' \code{\link{blackcsq3}} for Q_3 formula for futures.\\
#' \code{\link{blackcsprice}} for Corrado-Su option pricing formula for futures.
#'
#' @examples
#' blackcsq4(10, .3, 1.5, 0:100/100)
#' future <- 25
#' sigma <- 0.2
#' time <- 0.5
#' strike <- 26
#' rate <- 0.12
#' d1 <- blackd1(future, strike, time, rate, sigma)
#' blackcsq4(future, sigma, time, d1)
#'
#' @export
blackcsq4 <- function(future, sigma, time, d1){
q4 <- (d1 ^ 2 - 1 - 3 * sigma * sqrt(time) * (d1 - sigma * sqrt(time))) * dnorm(d1)
q4 <- q4 + (sigma * sigma * sigma * time ^ (1.5) * pnorm(d1))
q4 * future * sigma * sqrt(time) / 24
}
#' Computes call option price with Corrado-Su model for futures.
#'
#' Implementation of Corrado-Su formula for option pricing for futures.
#'
#' @param type 'call' for call option, any other value for put.
#' @param future current future price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param sigma volatility of the future price.
#' @param mu3 skewness.
#' @param mu4 kurtosis.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with option prices.
#'
#' @seealso
#' \code{\link{blackcsq3}} for Q_3 formula from Corrado-Su model for futures.\\
#' \code{\link{blackcsq4}} for Q_4 formula from Corrado-Su model for futures.
#'
#' @examples
#' spot <- 25
#' sigma <- 0.2
#' time <- 0.5
#' strike <- 26
#' rate <- 0.12
#' blackcsprice('call', spot, strike, time, rate, sigma, 0, 3) ==
#' blackprice('call', spot, strike, time, rate, sigma)
#' blackcsprice('call', spot, strike, time, rate, sigma, 0:100/100, 3)
#' blackcsprice('call', spot, strike, time, rate, sigma, 0.5, 0:1000/100)
#'
#' @export
blackcsprice <- function(type, future, strike, time, rate, sigma, mu3, mu4){
c_bs <- blackprice(type, future, strike, time, rate, sigma)
d1 <- (log(future / strike) + (sigma * sigma / 2) * time) / (sigma * sqrt(time))
q3 <- csq3(future, sigma, time, d1)
q4 <- csq4(future, sigma, time, d1)
callp <- c_bs + mu3 * q3 + (mu4 - 3) * q4
ifelse(tolower(type) == 'call', callp, callp - future + strike * exp(- rate * time))
}
#' Corrado-Su model fit for futures.
#'
#' Finds and returns Corrado-Su model parameters sigma, mu3 (skewness) and mu4 (kurtosis)
#' for a given set of options data.
#'
#' @param prices market prices for options
#' @param types 'call' for call option, any other value for put.
#' @param futures current future price.
#' @param strikes the strike price.
#' @param terms time to option expiration in years.
#' @param rates the risk-free interest rate.
#' @param initial_guess optional initial guess for optimization
#'
#' @return
#' List with model parameters sigma, mu3 and mu4
#'
#' @export
blackcs_fit <- function(prices, types, futures, strikes, terms, rates, initial_guess = NULL) {
lo <- c(0, -1, 1)
hi <- c(10, 1, 5)
if (is.null(initial_guess)) {
initial_guess <- c(CS.SIGMA, CS.SKEW, CS.KURT)
fit <- nloptr::nloptr(x0 = initial_guess,
eval_f = blackcs_minsq(prices, types, futures, strikes, terms, rates),
lb = lo,
ub = hi,
opts = list(algorithm = "NLOPT_GN_ISRES",
print_level = 0,
xtol_rel = 1e-4,
maxeval = 1000)
)
initial_guess <- fit$solution
}
fit <- nloptr::nloptr(x0 = initial_guess,
eval_f = blackcs_minsq(prices, types, futures, strikes, terms, rates),
lb = lo,
ub = hi,
opts = list(algorithm = "NLOPT_LN_COBYLA",
print_level = 0,
xtol_rel = 1e-4,
maxeval = 1000)
)
curr_param <- fit$solution
list(sigma = curr_param[1], mu3 = curr_param[2], mu4 = curr_param[3])
}
blackcs_minsq <- function(prices, types, futures, strikes, terms, rates) {
function(x) {
sigma <- x[1]
mu3 <- x[2]
mu4 <- x[3]
csp <- blackcsprice(types, futures, strikes, terms, rates, sigma, mu3, mu4)
if (any(is.nan(csp))) csp <- Inf
if (any(csp < 0)) csp[csp < 0] <- 0
sum((prices - csp)^2)
}
}
# used to estimate de implied volatility
# TODO: adaptar para corrado-su
cs.sigma.obj.func <- function(sigma, fut, strike, rate, time, mu3, mu4, price){
csprice(fut, strike, sigma, rate, time, mu3, mu4) - price
}
# TODO: adaptar para corrado-su -- por enquanto aplica apenas para calls
calc.sigma.cs <- function(batdf){
n <- length(batdf$Date)
sigma <- vector(mode='numeric', length=n)
for(i in 1:length(batdf$Date)){
root <- try(uniroot(f=cs.sigma.obj.func, interval=c(-0.1,1), tol=1e-15,
fut=batdf$Fut[i], strike=batdf$Strike[i],
rate=batdf$Rate[i], time=batdf$TBD[i],
mu3=batdf$mu3[i], mu4=batdf$mu4[i],
price=batdf$Opt.Price[i]))
if(inherits(root, "try-error")){
print(root)
print(paste("Inputs:", batdf$Fut[i], strike=batdf$Strike[i],
batdf$Rate[i], batdf$TBD[i], batdf$mu3[i], batdf$mu4[i],
batdf$Opt.Price[i], sep=' '))
print(paste('i=',i))
root <- list(root=-1)
}
sigma[i] <- root$root
}
return(sigma)
}
# Q14 56256
# V14 57151
# Z14 58056
# G15 58897
# M15 60635
# Z15 63594
#' Black Corrado-Su std. Deltas
#'
#' Finds and returns a dataframe with the strikes and implied
#' volatilites for the provided set of deltas or for the
#' default deltas.
#'
#' @param type Option type, "call" for calls or anything else for puts
#' @param future underlying future price at a given date
#' @param time time in years from a given date to a maturity date
#' @param rate interest free rate
#' @param sigma volatility of the stock price
#' @param mu3 skewness
#' @param mu4 kurtosis
#' @param deltas deltas for which strikes and vols should be calculated, 1-99. Defaults to usual
#' delta values.
#'
#' @return
#' Dataframe with columns delta, strike and impvol
#' @export
blackcs_deltas <- function(type, future, time, rate, sigma, mu3, mu4, deltas = DELTAS) {
eps <- 1e-8
strk_cs <- sapply(deltas/100, function(delta) {
uniroot(f=blackcs_strikes, interval=c(1, future*10), tol=eps,
delta=delta, type=type, future=future, time=time,
rate=rate, sigma=sigma, mu3=mu3, mu4=mu4)$root
})
impvol_cs <- sapply(strk_cs, function(strike) {
blackimpvol(blackcsprice(type, future, strike, time, rate, sigma, mu3, mu4),
type, future, strike, time, rate)
})
data.frame(delta=deltas, strike=strk_cs, impvol=impvol_cs)
}
blackcs_strikes <- function(x, delta, type, future, time, rate, sigma, mu3, mu4 ) {
csp <- blackcsprice(type, future, x, time, rate, sigma, mu3, mu4)
if (csp < 0) return(999)
varx <- tryCatch( { blackimpvol(csp, type, future, x, time, rate) },
error = function(e) {
1e-8
})
x - blackstrike(future, delta, varx, time)
}
#' Corrado-Su std. Deltas
#'
#' Finds and returns a dataframe with the strikes and implied
#' volatilites for the provided set of deltas or for the
#' default deltas.
#'
#' @param type Option type, "call" for calls or anything else for puts
#' @param spot underlying spot price at a given date
#' @param time time in years from a given date to a maturity date
#' @param rate interest free rate
#' @param yield dividend yield
#' @param sigma volatility of the stock price
#' @param mu3 skewness
#' @param mu4 kurtosis
#' @param deltas deltas for which strikes and vols should be calculated, 1-99. Defaults to usual
#' delta values.
#'
#' @return
#' Dataframe with columns delta, strike and impvol
#' @export
cs_deltas <- function(type, spot, time, rate, yield, sigma, mu3, mu4, deltas = DELTAS) {
ft <- spot * exp((rate - yield) * time)
eps <- 1e-8
strk_cs <- sapply(deltas/100, function(delta) {
uniroot(f=cs_strikes, interval=c(1, spot*10), tol=eps,
delta=delta, type=type, spot=spot, time=time,
rate=rate, yield=yield, sigma=sigma, mu3=mu3, mu4=mu4, ft=ft)$root
})
impvol_cs <-blackimpvol(csprice(type, spot, strk_cs, time, rate, yield, sigma, mu3, mu4),
type, ft, strk_cs, time, rate)
data.frame(delta=deltas, strike=strk_cs, impvol=impvol_cs)
}
cs_strikes <- function(x, delta, type, spot, strike, time, rate, yield, sigma, mu3, mu4, ft) {
csp <- csprice(type, spot, x, time, rate, yield, sigma, mu3, mu4)
if (csp < 0) return(999)
varx <- blackimpvol(csp, type, ft, x, time, rate)
x - blackstrike(ft, delta, varx, time)
}
var.garch <- function(rets, omega, alpha, beta) {
mu <- mean(rets)
e2 <- (rets - mu)^2
e2t <- omega + alpha*e2[c(-1,-length(e2))]
s2 <- stats::filter(e2t, beta, "recursive", init=mean(e2))
c(mean(e2), s2)
}
# source for student and ged formulas: https://hal.inria.fr/file/index/docid/270719/filename/B08027.pdf
.garch_likelihood <- function(rets, garchEnv, dist, scale=c(1, 1, 1)) {
.rets <- rets
.garchEnv <- garchEnv
llh_function <- switch(dist,
normal = {
function(...) {
x <- c(...)
omega <- x[1]/scale[1] # 1e-5
alpha <- x[2]/scale[2] # 1e-1
beta <- x[3]/scale[3]
garch <- var.garch(.rets, omega, alpha, beta)
# likelihood = 0.5 sum(-log(2 pi) - log(garch) - rets^2/garch)
# the implementation below is a simplification
if (any(garch < 0)) {
.llh <- get("llh", env=.garchEnv)
llh <- .llh + 0.1*(abs(.llh))
} else {
llh <- 0.5*sum(log(garch) + .rets[-1]^2/garch + log(2*pi))
}
assign("llh", llh, env=.garchEnv)
llh
}
},
student = {
function(...) {
x <- c(...)
omega <- x[1]/scale[1] # 1e-5
alpha <- x[2]/scale[2] # 1e-1
beta <- x[3]/scale[3]
v <- x[4]
n <- length(.rets)
garch <- var.garch(.rets, omega, alpha, beta)
s1 <- n*(lgamma((v + 1)/2) - lgamma(v/2) - (1/2)*log(pi*(v - 2)))
ll <- 1 + ((.rets[-1])^2/(garch * (v - 2)))
if (is.na(v) || any(garch < 0) || any(ll < 0)) {
.llh <- get("llh", env=.garchEnv)
llh <- .llh + 0.1*(abs(.llh))
} else {
llh <- - (s1 - 0.5*(sum(log(garch) + (v + 1)*log(ll))))
}
assign("llh", llh, env=.garchEnv)
llh
}
},
ged = {
function(...) {
x <- c(...)
omega <- x[1]/scale[1] # 1e-5
alpha <- x[2]/scale[2] # 1e-1
beta <- x[3]/scale[3]
v <- x[4]
n <- length(.rets)
yv <- sqrt((2^(-2/v)) * gamma(1/v) / gamma(3/v))
garch <- var.garch(.rets, omega, alpha, beta)
s1 <- n*(log(v/yv) - (1 + 1/v)*log(2) - lgamma(1/v))
if (is.na(v) || any(garch < 0)) {
.llh <- get("llh", env=.garchEnv)
llh <- .llh + 0.1*(abs(.llh))
} else {
llh <- - (s1 - 0.5*(sum(log(garch) + ((sqrt(garch))^(-v)) * ((abs(.rets[-1]/yv))^(v)))))
}
assign("llh", llh, env=.garchEnv)
llh
}
})
llh_function
}
garch_fit_1 <- function(rets, garchEnv, dist, bind_params) {
tiny <- 1e-3
alpha <- 0.1
beta <- 0.8
omega <- (1 - alpha - beta)*var(rets)
mu <- mean(rets)
v <- 1
x0 <- c(omega, alpha, beta)
# scaling and optimization bounds
small <- 1e-12
lo <- c(small, small, small)
hi <- c(2*var(rets), 1-small, 1-small)
# setting up likelihood function
garch_likelihood <- .garch_likelihood(rets, garchEnv, dist)
# optimization
eval_g0 <- function(x) {
c(-x[1], -x[2], -x[3], sum(x[2:3]) - (1 - tiny))
}
# student needs the degrees of freedom to be optimized too
if (dist == "student") {
v <- 3
x0 <- c(x0, v)
lo <- c(lo, 3)
hi <- c(hi, 10)
eval_g0 <- function(x) {
c(-x[1], -x[2], -x[3], sum(x[2:3]) - (1 - tiny), -x[4])
}
} else if (dist == "ged") { # ged needs nu to be estimated
x0 <- c(x0, v)
lo <- c(lo, 1)
hi <- c(hi, 10)
eval_g0 <- function(x) {
c(-x[1], -x[2], -x[3], sum(x[2:3]) - (1 - tiny), -x[4])
}
}
# optimization
if (bind_params) {
res <- nloptr::nloptr(x0=x0,
eval_f=garch_likelihood,
lb=lo,
ub=hi,
eval_g_ineq = eval_g0,
opts=list(algorithm="NLOPT_LN_COBYLA",
print_level = 0,
xtol_rel=1.0e-12, maxeval=2000))
} else {
res <- nloptr::nloptr(x0=x0,
eval_f=garch_likelihood,
lb=lo,
ub=hi,
opts=list(algorithm="NLOPT_LN_COBYLA",
print_level = 0,
xtol_rel=1.0e-12, maxeval=2000))
}
par <- res$solution
parms <- list(omega=par[1], alpha=par[2], beta=par[3])
var.g <- var.garch(rets, parms$omega, parms$alpha, parms$beta)
current_var <- sum(unlist(parms) * c(1, tail(rets, 1)^2, tail(var.g, 1)))
fit <- MyQuantTK::as.garchparms(c(parms, sigma_t=MyQuantTK::annualize(MyQuantTK::as.volatility(MyQuantTK::as.variance(current_var)))))
fit$shape <- par[4]
fit
}
.garch_likelihood_2 <- function(rets, garchEnv) {
.rets <- rets
.var <- as.numeric(var(rets))
.garchEnv <- garchEnv
function(...) {
x <- c(...)
omega <- .var*(1 - sum(x))
alpha <- x[1]
beta <- x[2]
garch <- var.garch(.rets, omega, alpha, beta)
# likelihood = 0.5 sum(-log(2 pi) - log(garch) - rets^2/garch)
# the implementation below is a simplification
if (any(garch < 0)) {
.llh <- get("llh", env=.garchEnv)
llh <- .llh + 0.1*(abs(.llh))
} else {
llh <- 0.5*sum(log(garch) + .rets[-1]^2/garch + log(2*pi))
}
assign("llh", llh, env=.garchEnv)
llh
}
}
garch_fit_2 <- function(rets, garchEnv, bind_params) {
tiny <- 1e-3
alpha <- 0.2
beta <- 0.7
mu <- mean(rets)
x0 <- c(alpha, beta)
# scaling and optimization bounds
small <- 1e-12
lo <- c(small, small)
hi <- c(1-small, 1-small)
# setting up likelihood function
garch_likelihood <- .garch_likelihood_2(rets, garchEnv)
# optimization
eval_g0 <- function(x) {
c(-x[1], -x[2], sum(x[1:2]) - (1 - tiny))
}
if (bind_params) {
res <- nloptr::nloptr(x0=x0,
eval_f=garch_likelihood,
lb=lo,
ub=hi,
eval_g_ineq = eval_g0,
opts=list(algorithm="NLOPT_LN_COBYLA",
print_level = 0,
xtol_rel=1.0e-12, maxeval=2000))
} else {
res <- nloptr::nloptr(x0=x0,
eval_f=garch_likelihood,
lb=lo,
ub=hi,
opts=list(algorithm="NLOPT_LN_COBYLA",
print_level = 0,
xtol_rel=1.0e-12, maxeval=2000))
}
par <- res$solution
parms <- list(omega=as.numeric(var(rets))*(1 - sum(par)), alpha=par[1], beta=par[2])
var.g <- var.garch(rets, parms$omega, parms$alpha, parms$beta)
current_var <- sum(unlist(parms) * c(1, tail(rets, 1)^2, tail(var.g, 1)))
fit <- MyQuantTK::as.garchparms(c(parms, sigma_t=MyQuantTK::annualize(MyQuantTK::as.volatility(MyQuantTK::as.variance(current_var)))))
fit$shape <- NA
fit
}
#' GARCH time series fitting
#'
#' Estimate parameters for a GARCH model, using Quasi-Maximum Likelihood Estimation
#' for optimization.
#'
#' @param rets log returns to be fitted
#' @param dist the conditional distribution to be used for estimation
#' @param bind_params if TRUE the fit will be bound by the alpha + beta < 1 restriction
#' @param force_asgard if TRUE will use asgard model instead of fair-garch
#'
#' @return
#' \code{"fit"} a list with the estimated omega, alpha and beta parameters,
#' and sigma, LTV, model, distribution and covariance matrix.
#'
#' @export
garch_fit <- function(rets, dist=c("normal", "student", "ged"), bind_params = TRUE, force_asgard = FALSE) {
dist <- match.arg(dist)
garchEnv <- new.env(hash = TRUE)
assign("llh", 1.0e99, env=garchEnv) # Keep the calculated likelihood function values to handle NaN
fit <- garch_fit_1(rets, garchEnv, dist, bind_params)
LTV <- MyQuantTK::ltv(fit)
model <- 'fair-garch'
if ((LTV < 0 || as.numeric(MyQuantTK::annualize(MyQuantTK::as.volatility(LTV))) > 500 ) || force_asgard) {
fit <- garch_fit_2(rets, garchEnv, bind_params)
model <- 'asgard'
}
LTV <- ltv(fit)
LTV_annual <- as.numeric(MyQuantTK::annualize(MyQuantTK::as.volatility(LTV)))
LTV <- as.numeric(LTV)
SIGMA_T <- as.numeric(fit$sigma_t)
fit$cvar <- cvar(rets, dist, fit$omega, fit$alpha, fit$beta, fit$shape)
ci <- conf_interval(fit$omega, fit$alpha, fit$beta, fit$cvar, LTV, LTV_annual)
list(omega = fit$omega, alpha = fit$alpha, beta = fit$beta, shape = fit$shape,
sigma_t = SIGMA_T, LTV = LTV, LTV_annual = LTV_annual,
model = model, dist = dist, cvar = fit$cvar, ci = ci)
}
cvar <- function(rets, dist, omega, alpha, beta, shape){
garchEnv <- new.env(hash = TRUE)
assign("llh", 1.0e99, env=garchEnv) # Keep the calculated likelihood function values to handle NaN
par <- c(omega, alpha, beta)
if(!is.na(shape))
par <- c(par, shape)
H <- hessian(.garch_likelihood(rets, garchEnv, dist), par,
method.args = list(eps = 1e-12))
H <- 0.5 * (H + t(H))
solve(H, tol = 1e-17)
}
conf_interval <- function(omega, alpha, beta, cvar, LTV, LTV_annual) {
dv_dw <- 1/(1 - alpha - beta)
dv_da <- omega/((1 - alpha - beta))^2
dv_db <- dv_da
var_vl <- (dv_dw ^ 2) * cvar[1,1] + (dv_da ^ 2) * cvar[2,2] +
(dv_db ^ 2) * cvar[3,3] + 2 * dv_dw * dv_da * cvar[1,2] +
2 * dv_dw * dv_db * cvar[1,3] + 2 * dv_da * dv_db * cvar[2,3]
.se_lt_var <- suppressWarnings(sqrt(var_vl))
volp <- MyQuantTK::as.volatility(LTV + .se_lt_var)[[1]]
volm <- MyQuantTK::as.volatility(LTV - .se_lt_var)[[1]]
if (!is.nan(volm) && !is.null(volm) && !is.na(volm) && volm < 0) volm <- 0
icmax <- sqrt(volp * 252)*100
icmin <- sqrt(volm * 252)*100
interv <- (icmax - icmin)/2
lo <- LTV_annual - interv
hi <- LTV_annual + interv
data.frame(Annualized_Volatility = LTV_annual,
CI_low = lo,
CI_high = hi,
In_CI = (LTV_annual <= hi && LTV_annual >= lo),
se_lt_var = .se_lt_var)
}
#' garchparms Class Creator
#'
#' Create list of class \code{"garchparms"}.
#'
#' @param omega Generalized AutoRegressive Conditional Heteroskedasticity
#' process paramater.
#' @param alpha GARCH process paramater.
#' @param beta GARCH process paramater.
#' @param sigma_t GARCH process paramater.
#'
#' @return
#' List with class \code{"garchparms"} from given arguments.
#'
#' @seealso
#' \code{\link{as.garchparms}}.
#'
#' @examples
#' garchparms(1, 1, 1, 1)
#' class(garchparms(1, 1, 1, 1))
#' typeof(garchparms(1, 1, 1, 1))
#'
#' @export
garchparms <- function(omega, alpha, beta, sigma_t) {
as.garchparms(list(
omega=omega,
alpha=alpha,
beta=beta,
sigma_t=sigma_t
))
}
#' garchparms Class Coercing
#'
#' Generic function for coercing to \code{"garchparms"}. There is no default
#' method, so only classes \code{"data.frame"}, \code{"numeric"}, \code{"list"}
#' are supported.
#'
#' @param x an object whose class will determine the method to be dispatched. Must
#' have the values for omega, alpha, beta, sigma_t.
#' @param ... further arguments to be passed to the next method.
#'
#' @return
#' List with class \code{garchparms} from given argument. When x has
#' \code{"data.frame"} class, idx = FALSE or equivalent will pass numeric(0) values
#' and numbers in idx will produce NA.
#'
#' @examples
#' class(garchparms(1, 1, 1, 1)) == class(list(omega=1, alpha=1, beta=1, sigma_t=1))
#' x = data.frame(omega=1, alpha=1, beta=1, sigma_t=1)
#' as.garchparms(x)
#' as.garchparms(x, 0)
#' as.garchparms(x, 2.579238)
#'
#' @name as.gachparms
#' @export
as.garchparms <- function(x, ...) {
UseMethod('as.garchparms', x)
}
#' @rdname as.gachparms
#'
#' @param idx additional parameter.
#'
#' @export
as.garchparms.data.frame <- function(x, idx=1, ...) {
as.garchparms.list(unclass(x[idx,]))
}
#' @rdname as.gachparms
#' @export
as.garchparms.numeric <- function(x, ...) {
as.garchparms.list(as.list(x))
}
#' @rdname as.gachparms
#' @export
as.garchparms.list <- function(x, ...) {
if (is.null(names(x)))
stop('names not defined')
aux <- c('omega', 'alpha', 'beta', 'sigma_t') %in% tolower(names(x))
if (!all(aux))
stop('names mismatch')
x$sigma_t <- as.volatility(x$sigma_t, annualized=T, dib=252)
structure(x, class='garchparms')
}
#' GARCH Information
#'
#' Gives relevant information about a GARCH model, such as parameters and
#' volatility. There is no default method, so only classes \code{"garch"},
#' \code{"fGARCH"}, \code{"uGARCHfit"} are supported.
#'
#' Methods for printing classes \code{"info.garch"}, \code{"info.fGARCH"},
#' \code{"info.uGARCHfit"} are also included.
#'
#' @param x a GARCH model.
#' @param ... further arguments to be passed to the next method.
#' @param dib days in one year, day count convention.
#'
#' @return
#' List with parameters, daily long term variance, annual long term variance, fit
#' of class "info.class(x)". Printing info classes will show daily and
#' annualized variance and volatility, plus parameters.
#'
#' @name info
#' @export
info <- function(x, ...) UseMethod('info')
#' @rdname info
#' @export
info.garch <- function(x, ..., dib=252) {
gamma <- ggamma(x)
lt_var <- coef(x)['omega'] / gamma
lt_vol <- sqrt(lt_var)
lt_vol_year <- lt_vol * sqrt(dib)
.vcov <- sqrt(diag(vcov(x)))
names(.vcov) <- c('omega', 'alpha', 'beta')
parms <- c(coef(x)['a0'], gamma, coef(x)['a1'], coef(x)['b1'])
names(parms) <- c('omega', 'gamma', 'alpha', 'beta')
.se_lt_var <- sqrt(sum(diag(vcov(x)) * c(1, lt_var ^ 2, lt_var ^ 2)) / (gamma ^ 2))
longterm <- c(lt_var, lt_vol, .se_lt_var)
names(longterm) <- c('variance', 'volatility', 'se')
longterm.annual <- c(lt_vol_year ^ 2, lt_vol_year)
names(longterm.annual) <- c('variance', 'volatility')
ans <- list(garch.parms=parms, longterm=longterm,
longterm.annual=longterm.annual,
fit=x, se=.vcov)
class(ans) <- 'info.garch'
ans
}
#' @rdname info
#' @export
print.info.garch <- function(x, ...) {
a <- x
class(a) <- NULL
cat('\nLong Term Info (%):', '\n')
print.default(a$longterm*rep(100, length(a$longterm)), digits=3,
print.gap=4)
cat('\nAnnualized Volatility (%):', '\n')
print.default(a$longterm.annual*rep(1, length(a$longterm.annual)),
digits=3, print.gap=4)
cat('\nGARCH parameters:', '\n')
print.default(a$garch.parms, digits=3, print.gap=4)
invisible(x)
}
#' @rdname info
#' @export
info.fGARCH <- function(x, ..., dib=252) {
# long term info
gamma <- ggamma(x)
lt_var <- ltv(x)
lt_vol <- as.volatility(lt_var)
lt_vol_year <- annualize(lt_vol, dib)
# se info
.vcov <- sqrt(diag(x@fit$cvar))
names(.vcov) <- c('omega', 'alpha', 'beta')
# parameters
parms <- c(coef(x)['omega'], gamma, coef(x)['alpha1'], coef(x)['beta1'])
names(parms) <- c('omega', 'gamma', 'alpha', 'beta')
# long term variance se
.se_lt_var <- sqrt(sum(diag(x@fit$cvar) * c(1, lt_var ^ 2, lt_var ^ 2)) / (gamma ^ 2))
longterm <- c(lt_var, lt_vol, .se_lt_var)
names(longterm) <- c('variance', 'volatility', 'se')
# annualized info
longterm.annual <- c(lt_vol_year ^ 2, lt_vol_year)
names(longterm.annual) <- c('variance', 'volatility')
# info structure
ans <- list(garch.parms=parms, longterm=longterm,
longterm.annual=longterm.annual,
fit=x, se=.vcov)
class(ans) <- 'info.garch'
ans
}
#' @rdname info
#' @export
print.info.fGARCH <- function(x, ...) {
print.info.garch(x, ...)
}
#' @rdname info
#' @export
info.uGARCHfit <- function(x, ..., dib=252) {
# long term info
gamma <- ggamma(x)
lt_var <- ltv(x)
lt_vol <- as.volatility(lt_var)
lt_vol_year <- annualize(lt_vol, dib)
# se info
.vcov <- sqrt(diag(x@fit$cvar))
names(.vcov) <- c('omega', 'alpha', 'beta')
# parameters
parms <- c(coef(x)['omega'], gamma, coef(x)['alpha1'], coef(x)['beta1'])
names(parms) <- c('omega', 'gamma', 'alpha', 'beta')
# long term variance se
.se_lt_var <- sqrt(sum(diag(x@fit$cvar) * c(1, lt_var ^ 2, lt_var ^ 2)) / (gamma ^ 2))
longterm <- c(lt_var, lt_vol, .se_lt_var)
names(longterm) <- c('variance', 'volatility', 'se')
# annualized info
longterm.annual <- c(lt_vol_year ^ 2, lt_vol_year)
names(longterm.annual) <- c('variance', 'volatility')
# info structure
ans <- list(garch.parms=parms, longterm=longterm,
longterm.annual=longterm.annual,
fit=x, se=.vcov)
class(ans) <- 'info.garch'
ans
}
#' @rdname info
#' @export
print.info.uGARCHfit <- function(x, ...) {
print.info.garch(x, ...)
}
#' GARCH Model Long Term Variance
#'
#' Calculates long term variance for GARCH models. There is no default
#' method, so only classes \code{"garch"}, \code{"fGARCH"}, \code{"uGARCHfit"},
#' \code{"garchparms"} are supported.
#'
#' @param x object from where omega, alpha and beta parameters will be extracted.
#'
#' @return
#' A structure with daily long term variance and class \code{"variance"}.
#'
#' @seealso
#' .
#'
#' @examples
#' ltv(garchparms(1, 2:4/20, .7, 1))
#'
#' @name ltv
#'
#' @export
ltv <- function(x) UseMethod('ltv')
#' @rdname ltv
#' @export
ltv.fGARCH <- function(x) as.variance(coef(x)['omega'] / ggamma(x))
#' @rdname ltv
#' @export
ltv.garch <- function(x) as.variance(coef(x)['a0'] / ggamma(x))
#' @rdname ltv
#' @export
ltv.garchparms <- function(x) as.variance(x$omega / ggamma(x))
#' @rdname ltv
#' @export
ltv.uGARCHfit <- function(x) as.variance(coef(x)['omega'] / ggamma(x))
#' GARCH Model Long Term Variance Decay
#'
#' Calculates long term variance decay for given alpha and beta. There is
#' no default method, so only classes \code{"fGARCH"}, \code{"garch"},
#' \code{"garchparms"}, \code{"uGARCHfit"} are supported.
#'
#' @param x object from where alpha and beta parameters will be extracted.
#'
#' @return
#' Numeric with long term variance decay for given parameters.
#'
#' @examples
#' decay(garchparms(1, 0:10/10, 1, 1))
#'
#' @name decay
#' @export
decay <- function(x) UseMethod('decay')
#' @rdname decay
#' @export
decay.fGARCH <- function(x) {
a <- log(1 / (coef(x)['alpha1'] + coef(x)['beta1']))
as.numeric(a)
}
#' @rdname decay
#' @export
decay.garch <- function(x) {
a <- log(1 / (coef(x)['a1'] + coef(x)['b1']))
as.numeric(a)
}
#' @rdname decay
#' @export
decay.garchparms <- function(x) {
a <- log(1 / (x$alpha + x$beta))
as.numeric(a)
}
#' @rdname decay
#' @export
decay.uGARCHfit <- function(x) {
a <- log(1 / (coef(x)['alpha1'] + coef(x)['beta1']))
as.numeric(a)
}
#' GARCH Model Gamma
#'
#' Computes gamma parameter for GARCH models. There is no default
#' method, so only classes \code{"garch"}, \code{"fGARCH"}, \code{"uGARCHfit"},
#' \code{"garchparms"} are supported.
#'
#' @param x an object whose class will determine the method to be dispatched.
#'
#' @return
#' Numeric with gamma calculated for given parameters.
#'
#' @seealso
#' .
#'
#' @examples
#' ggamma(garchparms(1, 2:4/20, .7, 1))
#'
#' @name ggamma
#' @export
ggamma <- function(x) UseMethod('ggamma')
#' @rdname ggamma
#' @export
ggamma.garch <- function(x) 1 - coef(x)['a1'] - coef(x)['b1']
#' @rdname ggamma
#' @export
ggamma.fGARCH <- function(x) 1 - coef(x)['alpha1'] - coef(x)['beta1']
#' @rdname ggamma
#' @export
ggamma.garchparms <- function(x) 1 - x$alpha - x$beta
#' @rdname ggamma
#' @export
ggamma.uGARCHfit <- function(x) 1 - coef(x)['alpha1'] - coef(x)['beta1']
#' GARCh Model Variance Term Structure
#'
#' Computes variance term structure for GARCH models. There is a default method,
#' where parameters for formula may be passed directly, and methods for classes
#' \code{"fGARCH"}, \code{"garchparms"}, \code{"uGARCHfit"}, where parameters are
#' calcuated from the GARCH models.
#'
#' @param x a GARCH model, GARCH paramters or ignored for default method.
#' @param ... further arguments to be passed to the next method.
#' @param v0 for default method, sigma_t converted to daily variance.
#' @param a for default method, long term variance decay.
#' @param vl for default method, long term variance.
#'
#' @return
#' Term structure of variance as a function of time.
#'
#' @seealso
#' \code{\link{decay}}, \code{\link{ltv}}.
#'
#' @examples
#' vts(garchparms(1, 1, 1, 1))(1)
#' vts("ignored", .5, .5, .3)(1)
#'
#' @name vts
#' @export
vts <- function(x, ...) UseMethod('vts')
#' @rdname vts
#' @export
vts.default <- function(x, v0, a, vl, ...) {
function(t) {
q <- (1 - exp(- a * t)) / (a * t)
as.variance(vl + q * (v0 - vl))
}
}
#' @rdname vts
#' @export
vts.fGARCH <- function(x, v0=NULL, ...) {
v0 <- if (is.null(v0)) tail(fBasics::volatility(x, type='h'), 1) else v0
a <- decay(x)
vl <- ltv(x)
vts.default(x, v0=v0, a=a, vl=vl)
}
#' @rdname vts
#' @export
vts.garchparms <- function(x, v0=NULL, ...) {
v0 <- if (is.null(v0)) daily(as.variance(x$sigma_t)) else v0
a <- decay(x)
vl <- ltv(x)
vts.default(x, v0=v0, a=a, vl=vl)
}
#' @rdname vts
#' @export
vts.uGARCHfit <- function(x, v0=NULL, ...) {
v0 <- if (is.null(v0)) daily(as.variance(x$sigma_t)) else v0
a <- decay(x)
vl <- ltv(x)
vts.default(x, v0=v0, a=a, vl=vl)
}
#' vts_update Volatility Term Structure Update
#'
#' Update volatility term structure given new return series data.
#'
#' @param garchdate current GARCH parameters estimation date.
#' @param refdate reference date for calculation.
#' @param returns return series to update vts.
#' @param omega omega GARCH parameter.
#' @param alpha alpha GARCH parameter.
#' @param beta beta GARCH parameter.
#' @param sigma_t sigma GARCH parameter.
#'
#' @return
#' Volatility Term Structure.
#'
#' @name vts_update
#' @export
vts_update <- function(garchdate, refdate, returns, omega, alpha, beta, sigma_t) {
if(garchdate < refdate) {
sigma_t <- sigma_t_update(returns, omega, alpha, beta, sigma_t)
sigma_t <- as.numeric(sigma_t)
}
vts(garchparms(omega, alpha, beta, sigma_t))
}
# sigma(t) ----
#' sigma_t Calculation for GARCH Models
#'
#' Calculates sigma_t for GARCH models. There is no default method, so only
#' classes \code{"fGARCH"}, \code{"uGARCHfit"}, are supported.
#'
#' @param x a GARCH model with class \code{"fGARCH"} or \code{"uGARCHfit"}.
#' @param ... further arguments to be passed to the next method.
#'
#' @return
#' sigma_t for GARCH model in annualized volatility.
#'
#' @name sigma_t
#' @export
sigma_t <- function(x, ...) UseMethod('sigma_t')
#' @rdname sigma_t
#' @export
sigma_t.fGARCH <- function(x, ...) {
last_ret <- tail(x@data, 1)
last_var <- tail(x@sigma.t, 1) ^ 2
current_var <- sum(coef(x) * c(1, last_ret ^ 2, last_var))
annualize(as.volatility(as.variance(current_var)))
}
# #' @rdname sigma_t
# #' @export
# sigma_t.uGARCHfit <- function(fit) {
# current_var <- rugarch::sigma(rugarch::ugarchforecast(fit, n.ahead = 1)) ^ 2
# annualize(as.volatility(as.variance(current_var)))
# }
#' sigma_t_update Update GARCH paramater Sigma
#'
#' Given an additional return series, update sigma GARCH parameter, keeping the
#' others.
#'
#' @param returns return series to update sigma.
#' @param omega omega Garch parameter.
#' @param alpha alpha Garch parameter.
#' @param beta beta Garch parameter.
#' @param sigma_t sigma Garch parameter previously estimated.
#'
#' @return
#' Annualized sigma_t.
#'
#' @name sigma_t_update
#' @export
sigma_t_update <- function(returns, omega, alpha, beta, sigma_t) {
sig <- daily(as.volatility(sigma_t, annualized = T, dib = 252)) ^ 2
x <- omega + alpha * returns ^ 2
if(length(x) == 1)
sig <- x + beta * sig
else
sig <- filter(x, beta, method = "recursive", init = sig)
sig <- sig[length(sig)]
annualize(as.volatility(as.variance(sig)))
}
# volatility ----
#' volatility_update Volatility Estimative Update
#'
#' Estimate new volatility value within initial date with known volatility
#' an desired final date.
#'
#' @param ticker underlying ticker identification.
#' @param refdate reference date for calculation.
#' @param matdate desired final date to witch update volatility.
#' @param garchdate current GARCH parameters estimation date.
#' @param omega omega GARCH parameter.
#' @param alpha alpha GARCH parameter.
#' @param beta beta GARCH parameter.
#' @param sigma_t sigma GARCH parameter.
#'
#' @return
#' Annualized volatility.
#'
#' @name volatility_update
#' @export
volatility_update <- function(ticker, refdate, matdate, garchdate, omega, alpha, beta, sigma_t) {
vts <- vts_update(ticker, garchdate, refdate, omega, alpha, beta, sigma_t)
bizdays <- bizdays::bizdays(refdate, matdate)
volatility <- if (! is.null(vts)) vts(bizdays) else NA
as.numeric(annualize(as.volatility(volatility))) / 100
}
#' instrate Continuous Compounded Interest rate interpolation
#'
#' Given an interest rate curve, find the continuous compounded interest
#' rate for an instant in time.
#'
#' @param curve yield curve.
#' @param refdate reference date for calculation.
#' @param matdate desired maturity for interest rate compounding.
#'
#' @return
#' Continuously compounded interest rate.
#'
#' @name instrate
#' @export
instrate <- function(curve, refdate, matdate) {
bizdays <- bizdays::bizdays(refdate, matdate)
rate <- curve[bizdays]
as.numeric(log(1 + rate / 100))
}
# #' garchparms_info Garch fitting parameters summary.
# #'
# #' Summarize parameters omega, alpha, beta, sigma and long-term volatility
# #' from input serie.
# #'
# #' @param .fit valid GARCH fitting result.
# #' @param series input series for fitted data.
# #'
# #' @return
# #' Dataframe with parameters omega, alpha, beta, sigma and ltv in this order.
# #'
# #' @name garchparms_info
# #' @export
# garchparms_info <- function(.fit, series) {
# if(is.null(series)) return(NULL)
# series <- series[!is.na(series)]
#
# fit <- tryCatch({ .fit(series) }, error=function(e) { NULL })
# if (is.null(fit)) {
# warning(ticker, ': nao foi possivel realizar o fit!')
# return(NULL)
# }
#
# fit.info <- info(fit)
# lt_vol <- as.numeric(annualize(as.volatility(ltv(fit))))
# if(is.infinite(lt_vol) | is.na(lt_vol)) warning(ticker, ': LTV nao converge!')
#
# l <- list(omega=fit.info$garch.parms['omega'],
# alpha=fit.info$garch.parms['alpha'],
# beta=fit.info$garch.parms['beta'],
# 'sigma_t'=as.numeric(sigma_t(fit)),
# 'LTV'=lt_vol)
# as.data.frame(l)
# }
#
# #' @export
# garchparms_calc <- function(ticker, from=NULL, to, gpack=c("fGarch", "rugarch")) {
# gpack <- match.arg(gpack)
# .fit <- switch(gpack,
# fGarch = function(x) {
# fGarch::garchFit(~garch(1,1), x, cond.dist='QMLE',
# include.mean=FALSE, trace=FALSE)
# },
# rugarch = function(x) {
# gspec.ru <- rugarch::ugarchspec(
# variance.model = list(model="sGARCH", garchOrder=c(1, 1)),
# mean.model=list(armaOrder=c(0,0), include.mean=F),
# distribution.model="norm")
# rugarch::ugarchfit(gspec.ru, x)
# }
# )
#
# to <- as.Date(bizdays::adjust.previous(to))
# if(is.null(from)) from <- to - lubridate::dyears(3)
# from <- as.Date(bizdays::adjust.next(from))
#
# df <- data.frame(ticker, from, to, stringsAsFactors = F)
# garch <- apply(df, 1, function(x, .fit) {
# garchparms_info(.fit, x['ticker'], x['from'], x['to'])
# }, .fit=.fit)
#
# if(!is.null(garch)) garch <- do.call(rbind, garch)
# row.names(garch) <- NULL
# garch
# }
#' gmoa.pricing
#'
#' BVMF pricing models package.
#'
#' @name gmoa.pricing-package
#' @import bizdays
#' @import dplyr
#' @importFrom numDeriv hessian
#' @importFrom fBasics volatility
#' @importFrom reshape melt
NULL
#' Computes european call option price via Heston Model
#'
#' Heston formula implementation, pricing european call options.
#'
#' @param S spot price
#' @param k the strike price.
#' @param tau annualized period between option reference date and maturity date
#' @param r interest free rate
#' @param q dividend yield rate
#' @param lambda speed of the reversion mean of the volatility
#' @param eta volatility of the volatility
#' @param rho correlation between the Brownian motion of St and vt.
#' @param vbar long term mean of volatility
#' @param v0 initial variance
#'
#' @return
#' Vector with price of european options, according to
#' Heston model.
#'
#' @export
heston_price <- function(S, K, tau, r, q, lambda, eta, rho, vbar, v0) {
f <- S * exp((r - q) * tau)
x <- log(f/K)
integral <- heston_phi_transform(tau, x, lambda, eta, rho, vbar, v0)
f - (sqrt(K*f)/(2*pi)) * integral
}
#' Computes european call option price via Heston Model
#'
#' Heston formula implementation, pricing european call futures options.
#'
#' @param F future price
#' @param k the strike price.
#' @param tau annualized period between option reference date and maturity date
#' @param lambda speed of the reversion mean of the volatility
#' @param eta volatility of the volatility
#' @param rho correlation between the Brownian motion of St and vt.
#' @param vbar long term mean of volatility
#' @param v0 initial variance
#'
#' @return
#' Vector with price of european options, according to
#' Heston model.
#'
#' @export
heston_pricef <- function(f, K, tau, lambda, eta, rho, vbar, v0) {
x <- log(f/K)
integral <- heston_phi_transform(tau, x, lambda, eta, rho, vbar, v0)
f - (sqrt(K*f)/(2*pi)) * integral
}
#' Heston model fit
#'
#' Fits a heston model to the option data provided and returns
#' the model parameters.
#'
#' @param forwards forward prices of the underlying
#' @param strikes strike prices of the options
#' @param prices options closing prices for the reference dates
#' @param terms annualized periods between option reference dates and maturity dates
#' @param initial_guess vector with the initial guess for the parameters lambda, eta, rho, vbar and v0
#' @param uncertainties price values uncertainties in \% for weighting
#' @param do_imp_vol if TRUE will minimize implied vol error (market imp vols must be passed in prices param)
#' @param do_global if TRUE does global optimization to find starting point for local optimization
#'
#' @return
#' List with heston model parameters lambda, eta, rho, vbar and v0
#'
#' @export
heston_fitf <- function(forwards, strikes, prices, terms,
initial_guess = NULL, uncertainties = NULL,
do_imp_vol = FALSE, do_global = TRUE) {
min_function <- heston_minsqf(forwards, strikes, prices, terms, uncertainties, do_imp_vol)
heston_optimize(min_function, initial_guess, do_global)
}
heston_minsqf <- function(forwards, strikes, prices, terms, uncertainties = NULL, do_imp_vol = FALSE) {
function(...) {
x <- c(...)
lambda <- x[1]
eta <- x[2]
rho <- x[3]
vbar <- x[4]
v0 <- x[5]
heston_prices <- numeric()
# Integration doesn't work well for vectors
for(i in 1 : length(forwards)) {
heston_prices[i] <- heston_pricef(forwards[i], strikes[i], terms[i],
lambda, eta, rho, vbar, v0)
}
if (do_imp_vol) {
heston_prices <- tryCatch({
gmoa.pricing::bsmimpvol(heston_prices, "call", spots, strikes, terms, rates, yields)
}, error = function(e) { 0 })
}
if (is.null(uncertainties)) {
sum(((heston_prices - prices)/prices)^2)
} else {
sum(((((heston_prices - prices)/prices)) * 1/((uncertainties/100)))^2)
}
}
}
#' Heston model fit
#'
#' Fits a heston model to the stock option data provided and returns
#' the model parameters.
#'
#' @param spots forward prices of the underlying
#' @param strikes strike prices of the options
#' @param prices options closing prices for the reference dates
#' @param rates interest free rates
#' @param yields dividend yields
#' @param terms annualized periods between option reference dates and maturity dates
#' @param initial_guess vector with the initial guess for the parameters lambda, eta, rho, vbar and v0
#' @param uncertainties price values uncertainties in \% for weighting
#' @param do_imp_vol if TRUE will minimize implied vol error (market imp vols must be passed in prices param)
#' @param do_global if TRUE does global optimization to find starting point for local optimization
#'
#' @return
#' List with heston model parameters lambda, eta, rho, vbar and v0
#'
#' @export
heston_fit <- function(spots, strikes, prices, rates, yields, terms,
initial_guess = NULL, uncertainties = NULL,
do_imp_vol = FALSE, do_global = TRUE) {
min_function <- heston_minsq(spots, strikes, prices, rates, yields, terms, uncertainties, do_imp_vol)
heston_optimize(min_function, initial_guess, do_global)
}
heston_minsq <- function(spots, strikes, prices, rates, yields, terms, uncertainties = NULL, do_imp_vol = FALSE) {
function(...) {
x <- c(...)
lambda <- x[1]
eta <- x[2]
rho <- x[3]
vbar <- x[4]
v0 <- x[5]
heston_prices <- numeric()
# Integration doesn't work well for vectors
for(i in 1 : length(spots)) {
heston_prices[i] <- heston_price(spots[i], strikes[i], terms[i], rates[i],
yields[i], lambda, eta, rho, vbar, v0)
}
if (do_imp_vol) {
heston_prices <- tryCatch({
gmoa.pricing::bsmimpvol(heston_prices, "call", spots, strikes, terms, rates, yields)
}, error = function(e) { 0 })
}
if (is.null(uncertainties)) {
sum(((heston_prices - prices)/prices)^2)
} else {
sum(((((heston_prices - prices)/prices)) * 1/((uncertainties/100)))^2)
}
}
}
heston_phi <- function(k, tau, lambda, eta, rho, vbar, v0) {
b <- lambda + complex (imaginary = rho * eta) * k
d <- sqrt(b^2 + (eta^2) * (k) * (k + complex(imaginary = -1)))
g <- (b - d)/(b + d)
T_m <- (b - d)/(eta ^ 2)
t <- T_m * ( 1 - exp(-d * tau) )/( 1 - g * exp(-d * tau) )
W <- lambda * vbar * (tau * T_m - 2 * log(( 1 - g*exp(-d * tau) )/( 1 - g ))/(eta ^ 2))
exp(W + v0*t)
}
heston_phi_transform <- function(tau, x, lambda, eta, rho, vbar, v0) {
integrand <- function(k) {
result <- 2 * Re(exp(complex(imaginary = -k * x)) * heston_phi(complex(real = k, imaginary = 0.5), tau,
lambda, eta, rho, vbar, v0))/(k ^ 2 + 1/4)
if (!all(is.numeric(result))) browser()
if (any(is.nan(result))) browser()
if (any(is.na(result))) browser()
result
}
integrate(integrand, 0, 100, stop.on.error = F)$value
}
heston_optimize <- function(objective_function, initial_guess = NULL, do_global = FALSE) {
if (is.null(initial_guess)) {
x0 <- c(1e-6, 1e-6, -1, 1e-6, 1e-6)
} else {
x0 <- initial_guess
}
lo <- c(1e-12, 1e-12, -1, 1e-12, 1e-12)
hi <- c(20, 5, 0, 1, 1)
# Feller's condition
# 2*lambda*vbar - eta^2 > 0
# lambda, eta, rho, vbar, v0
eval_g0 <- function(x) {
c(x[2]*x[2] - (2*x[1]*x[4]))
}
if (do_global) {
gl_fit <- nloptr::nloptr(x0 = x0,
eval_f = objective_function,
lb = lo,
ub = hi,
eval_g_ineq = eval_g0,
opts = list(algorithm="NLOPT_GN_ISRES",
print_level = 0,
xtol_rel=1.0e-6, maxeval = 500))
x0 <- gl_fit$solution
}
fit <- nloptr::nloptr(x0 = x0,
eval_f = objective_function,
lb = lo,
ub = hi,
eval_g_ineq = eval_g0,
opts = list(algorithm="NLOPT_LN_COBYLA",
print_level = 0,
xtol_rel=1.0e-6, maxeval = 500))
params <- fit$solution
list(lambda = params[1], eta = params[2], rho = params[3],
vbar = params[4] , v0 = params[5], objective = fit$objective)
}
#' Heston model std. Deltas
#'
#' Finds and returns a dataframe with the strikes and implied
#' volatilites for the provided set of deltas or for the
#' default deltas.
#'
#' @param fit a Heston model fit as returned by heston_fit
#' @param S underlying spot price at a given date
#' @param tau time in years from a given date to a maturity date
#' @param r interest free rate
#' @param q dividend yield
#' @param future underlying future price, can be provided instead of spot and yield
#' @param deltas deltas for which strikes and vols should be calculated. Defaults to usual
#' delta values.
#'
#' @return
#' Dataframe with columns delta, strike and impvol
#'
#' @export
heston_deltas <- function(fit, S = NULL, tau, r, q = NULL, future = NULL, deltas = DELTAS) {
if (is.null(future)) {
ft <- S * exp((r - q) * tau)
} else {
ft <- future
}
if (is.null(ft)) stop("Future price or spot, rate and yield must be provided.")
eps <- 1e-8
# heston price não funciona vetorialmente devido à integral,
# então as chamadas precisam ficar em apply
strk_heston <- sapply(deltas/100, function(delta) {
uniroot(f=heston_strikes, interval=c(1, ft*10), tol=eps,
t=tau, delta=delta, ft=ft, fit=fit, r=r)$root
})
impvol_heston <- sapply(strk_heston, function(strike) {
blackimpvol(heston_pricef(ft, strike, tau, fit$lambda, fit$eta, fit$rho, fit$vbar, fit$v0), "call", ft, strike, tau, r)
})
data.frame(delta=deltas, strike=strk_heston, impvol=impvol_heston)
}
heston_strikes <- function(x, delta, ft, t, fit, r) {
hp <- heston_pricef(ft, x, t, fit$lambda, fit$eta, fit$rho, fit$vbar, fit$v0)
if (hp < 0) return(999)
varx <- tryCatch( { blackimpvol(hp, "call", ft, x, t, r) },
error = function(e) {
1e-8
})
x - blackstrike(ft, delta, varx, t)
}
#' Computes implied volatility inverting Black-Scholes-Merton formula.
#'
#' Uses bisection method to calculate implied volatility, given other parameters
#' for Black-Scholes-Merton formula for pricing options on stocks.
#'
#' @param option_prices stock option price.
#' @param type 'call' for call option, any other value for put.
#' @param spot current stock price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param yield the dividends that are expected to be paid.
#' @param tolerance the tolerance for the approximation.
#' @param maxiter maximum number of iterations for bisection method.
#' @param check function for checking the convergence of all bisections.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' An approximation for the implied volatility. If the algorithm converges before
#' \code{maxiter} is reached, Black-Scholes-Merton formula calculated with this
#' volatility should not differ from \code{option_prices} by more than
#' \code{tolerance}.
#'
#' @seealso
#' \code{\link{blackimpvol}} for options on futures, with Black-76 model.
#'
#' @examples
#' bsmimpvol(6, 'call', 50, 52, 1, 0.1, 0)
#' bsmimpvol(c(3, 4, 5, 6), 'call', 50, 52, 1, 0.1, 0)
#' bsmimpvol(1.5, 'put', 40, 38, 0.5, 0.15, 0, tolerance = 1e-8)
#'
#' @export
bsmimpvol <- function(option_prices, type, spot, strike, time, rate, yield,
tolerance=.Machine$double.eps, maxiter=100, check=any) {
f <- function(sigma, ...){
bsmprice(sigma=sigma, ...) - option_prices
}
res <- multiroot(func_=f, interval=c(1e-8, 10), tolerance=tolerance,
maxiter=maxiter, check=check, type=type, spot=spot,
strike=strike, rate=rate, time=time, yield=yield)
res$root
}
#' Computes implied volatility inverting Black-76 formula.
#'
#' Uses bisection method to calculate implied volatility, given other parameters
#' for Black-76 formula for pricing options on futures.
#'
#' @param option_prices future option price.
#' @param type 'call' for call option, any other value for put option.
#' @param future current future price.
#' @param strike the strike price.
#' @param time time to option expiration in years.
#' @param rate the risk-free interest rate.
#' @param tolerance the tolerance for the approximation.
#' @param maxiter maximum number of iterations for bisection method.
#' @param check function for checking the convergence of all bisections.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with approximations implied volatility. If the algorithm converges
#' before \code{maxiter} is reached, Black-76 formula calculated with this
#' volatility should not differ from \code{option_prices} by more than
#' \code{tolerance}.
#'
#' @seealso
#' \code{\link{bsmimpvol}} for options on stocks, with Black-Scholes-Merton model.
#'
#' @examples
#' blackimpvol(c(3.66765, 5.441491), c('call', 'put'), 50, 52, 1, 0.12)
#' blackimpvol(c(2, 3, 4, 5, 6, 7), 'put', 50, 52, 1, 0.12)
#' blackimpvol(2000:6000/1000, 'call', 50, 52, 1, 0.12)
#'
#' @export
blackimpvol <- function(option_prices, type, future, strike, time, rate,
tolerance=.Machine$double.eps, maxiter=100, check=any) {
f <- function(sigma, ...){
blackprice(sigma=sigma, ...) - option_prices
}
res <- multiroot(func_=f, interval=c(1e-8, 10), tolerance=tolerance,
maxiter=maxiter, check=check, type=type, future=future,
strike=strike, rate=rate, time=time)
res$root
}
#' Vectorial bisection method varying parameters from one function.
#'
#' Searches the interval for a root (i.e., zero) of the function func_ with
#' respect to its first argument. Other arguments may be vectorial, so
#' \code{multiroot} will return a vector with roots.
#'
#' \code{multiroot} was created aiming implied volatility calculation for
#' several options.
#'
#' @param func_ objective function for root finding.
#' @param interval a vector containing the end-points of the interval to be
#' searched for all the roots.
#' @param ... additional named or unnamed arguments to be passed to func_,
#' may be vectorial, and this is the difference between multiroot and uniroot.
#' @param tolerance the desired accuracy (convergence tolerance).
#' @param maxiter the maximum number of iterations.
#' @param check function for checking convergence in \code{\link{any}} or
#' \code{\link{all}} of the roots. Since it is bisection method, they are
#' the same.
#'
#' @return
#' Vector with roots.
#'
#' @seealso
#' \code{\link{uniroot}} for R built-in one dimensional root finding.
#'
#' @examples
#' f <- function(x, y) x - y
#' multiroot(f, c(-3, 10), 1:3)
#' g <- function(sigma, premium, type, spot, strike, time, rate, yield)
#' bsmprice(type, spot, strike, time, rate, yield, sigma) - premium
#' multiroot(g, c(1e-8, 10), 5, 'call', 30, 31, c(1, 1.1, 1.2, 1.3), 0.15, 0)
#' premium <- blackprice('put', 10, 10, 0.5, 0.12, 1:100/100)
#' h <- function(sigma, premium, type, spot, strike, time, rate)
#' blackprice(type, spot, strike, time, rate, sigma) - premium
#' multiroot(h, c(1e-8, 10), premium, 'put', 10, 10, 0.5, 0.12)
#'
#' @export
multiroot <- function(func_, interval, ...,
tolerance=.Machine$double.eps, maxiter=100, check=any) {
s_lower <- sign(func_(interval[1], ...))
s_upper <- sign(func_(interval[2], ...))
if ( any(s_upper*s_lower != - 1) )
stop('function evaluated at upper and lower limits must have the opposite sign')
func <- function(a, ...) func_(a, ...) * s_upper
initial_guess <- sum(interval) / 2
x <- initial_guess * rep(1, length(s_upper))
lower <- interval[1] * rep(1, length(s_upper))
upper <- interval[2] * rep(1, length(s_upper))
iter <- 0
err <- func(x, ...)
while (check(abs(err) > tolerance) && iter < maxiter) {
idx.err <- err < 0
idx.err <- ifelse(is.na(idx.err), FALSE, idx.err)
lower[idx.err] <- x[idx.err]
x[idx.err] <- (upper[idx.err] + x[idx.err]) / 2
idx.err <- ! idx.err
upper[idx.err] <- x[idx.err]
x[idx.err] <- (lower[idx.err] + x[idx.err]) / 2
err <- func(x, ...)
iter <- iter + 1
}
res <- list(root=x, iter=iter, err=err)
structure(res, class='multiroot')
}
#' Print Bisection Method Result.
#'
#' Specific method for printing the multiroot class, created for multiroot
#' function, in a pleasing form.
#'
#' The multiroot class is a list with slots root, iter, err.
#'
#' @param x a list.
#' @param ... additional arguments, will be ignored.
#'
#' @return
#' \code{invisible(x)}.
#'
#' @seealso
#' \code{\link{multiroot}} for vectorial bisection method.
#'
#' @examples
#' mt <- structure(list(root=1, iter=20, err=1e-8), class='multiroot')
#' print(mt)
#'
#' @export
print.multiroot <- function(x, ...) {
cat('root:', x$root, '\n')
cat('iter:', x$iter, '\n')
cat('err:', x$err, '\n')
invisible(x)
}
#' z auxiliary function for SABR(Stochastic Alpha Beta Rho) volatility model.
#'
#' SABR model formula is broken in auxiliary functions, and z is one
#' of them.
#'
#' @param fut current future price.
#' @param strike the strike price.
#' @param alpha parameter from model.
#' @param beta parameter from model.
#' @param nu parameter from model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with the z formula values.
#'
#' @seealso
#' \code{\link{x.func}} for x auxiliary function from SABR model.
#' \code{\link{sabrimpvol}} for volatility fit with SABR model.
#'
#' @examples
#' z.func(10, 9, 0.5, 0.5, 0.5)
#' z.func(50:150/10, 9, 0.5, 0.5, 0.5)
#' z.func(50, 45, 0.4325, 0.65846, 1.5)
#'
#' @export
z.func <- function(fut, strike, alpha, beta, nu){
return((nu / alpha) * ((fut * strike) ^ (0.5 * (1 - beta))) * log(fut / strike))
}
#' x auxiliary function for SABR(Stochastic Alpha Beta Rho) volatility model.
#'
#' SABR model formula is separated in auxiliary functions, and x is one
#' of them.
#'
#' @param z other auxiliary function for SABR.
#' @param rho parameter from model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with the x formula values.
#'
#' @seealso
#' \code{\link{z.func}} for z auxiliary function from SABR model.
#' \code{\link{sabrimpvol}} for volatility fit with SABR model.
#'
#' @examples
#' x.func(z.func(10, 9, 0.5, 0.5, 0.5), 0.5)
#' x.func(z.func(50, 45, 0.4325, 0.65846, 1.5), -0.15678)
#' x.func(-100:100/100, 0)
#'
#' @export
x.func <- function(z, rho){
fct <- 1 - 2 * rho * z + z ^ 2
fct <- ifelse(fct < 0, 0, fct)
fct <- sqrt(fct)
#den <- 1 - rho
log.arg <- (fct + z - rho) / (1 - rho)
log.arg <- ifelse(log.arg < 0, 0, log.arg)
#TODO: verificar caso em que 1 - rho == 0
return( log(log.arg) )
}
#' Volatility fit with SABR(Stochastic Alpha Beta Rho) model.
#'
#' SABR model for volatility market data fit implementation.
#'
#' @param spot current stock price.
#' @param rate the risk-free interest rate.
#' @param strike the strike price.
#' @param bizdays time to option expiration in years, based in working days.
#' @param alpha parameter from model.
#' @param beta parameter from model.
#' @param rho parameter from model.
#' @param nu parameter from model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with fit volatility, according to SABR model.
#'
#' @seealso
#' \code{\link{z.func}} for z auxiliary function from SABR model.
#' \code{\link{x.func}} for x auxiliary function from SABR model.
#'
#' @examples
#' sabrimpvol(5, 0.16, 4.5, 1, 0.5, 0.5, 0.5, 0.5)
#' sabrimpvol(5, 0.16, 4.5, 1, 0.5, 0:100/100, 0.5, 0.5)
#' sabrimpvol(5, .16, 4.5, 1, .5, .5, -100:100/100, .5)
#'
#' @export
sabrimpvol <- function(spot = NULL, rate = NULL, fut = NULL, strike, bizdays, alpha, beta, rho, nu){
if (is.null(fut)) {
fut <- spot * exp(rate * bizdays)
}
z <- z.func(fut, strike, alpha, beta, nu)
x <- x.func(z, rho)
dn1 <- (fut * strike) ^ (0.5 * (1 - beta))
dn2 <- 1 + (1 / 24) * ((1 - beta) ^ 2) * ((log(fut / strike)) ^ 2) +
(1 / 1920) * ((1 - beta) ^ 4) * (log(fut / strike)) ^ 4
term1 <- (alpha / (dn1 * dn2)) * (z / x)
fct <- (1 / 24) * ((1 - beta) ^ 2) * (alpha ^ 2) / ((fut * strike) ^ (1 - beta)) +
(1 / 4) * (rho * beta * nu * alpha) / ((fut * strike) ^ (0.5 * (1 - beta))) +
(1 / 24) * (2 - 3 * (rho ^ 2)) * nu ^ 2
term2 <- 1 + (fct * bizdays)
impvol <- term1 * term2
if (any(is.nan(impvol))) impvol[is.nan(impvol)] <- 0
return(impvol)
}
#' Numerical derivative of implied volatility formula from SABR method with
#' respect to alpha.
#'
#' Symmetric Difference Quotient method of numerical differentiation applied
#' in the volatility formula from SABR method, with respect to alpha.
#'
#' @param spot current stock price.
#' @param rate the risk-free interest rate.
#' @param strike the strike price.
#' @param bizdays time to option expiration in years, based in working days.
#' @param alpha parameter from model.
#' @param beta parameter from model.
#' @param rho parameter from model.
#' @param nu parameter from model.
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Vector with the derivatives.
#'
#' @seealso
#' \code{\link{sabrimpvol}} for SABR volatility formula.
#'
#' @examples
#' dsabrimpvol.dalpha(5, 0.16, 4.5, 1, 0.5, 0.5, 0.5, 0.5)
#' dsabrimpvol.dalpha(5, 0.16, 4.5, 1, 0.5, 0:100/100, 0.5, 0.5)
#' dsabrimpvol.dalpha(5, .16, 4.5, 1, .5, .5, -100:100/100, .5)
#'
#' @export
dsabrimpvol.dalpha <- function(spot, rate, strike, bizdays, alpha, beta, rho, nu){
d.alpha <- 1e-6
impvol.up <- sabrimpvol(spot, rate, strike, bizdays, alpha + d.alpha, beta, rho, nu)
impvol.down <- sabrimpvol(spot, rate, strike, bizdays, alpha - d.alpha, beta, rho, nu)
return((impvol.up - impvol.down) / (2 * d.alpha))
}
#' SABR model fit
#'
#' Finds and returns SABR model parameters alpha, beta, nu and rho for
#' a given set of strikes and sigmas.
#'
#' @param spot current stock price. If not provided will try to use the future parameter.
#' @param rate the risk-free interest rate. If not provided will try to use the future parameter.
#' @param future the future price. If not provided will calculate future from spot and rate.
#' @param tau time to option expiration in years, based in working days.
#' @param strikes options strikes
#' @param sigmas options sigmas
#' @param initial_guess named vector in order alpha, nu, rho, beta of initial parameters for optimization (optional)
#'
#' @return
#' List with fitted parameters alpha, beta, nu, rho
#'
#' @seealso
#' \code{\link{sabrimpvol}} for SABR volatility
#'
#' @export
sabr_fit <- function(future = NULL, strikes, tau, sigmas, spot = NULL, rate = NULL, initial_guess = NULL) {
if (is.null(future) && (is.null(spot) || is.null(rate))) stop("Either future price or spot and rate must be provided.")
if (is.null(future)) future <- spot * exp(rate * tau)
if (is.null(initial_guess)) initial_guess <- c(alpha=ALPHA, nu=NU, rho=RHO, beta=BETA)
fit <- nloptr::lbfgs(x0=initial_guess,
fn=sabr_optim(future, strikes, tau, sigmas),
lower=c(0.0001, 0.0001, -1, 0),
upper=c(20, Inf, 0.99, 1))
list(alpha = fit$par[1], nu = fit$par[2], rho = fit$par[3], beta = fit$par[4])
}
sabr_optim <- function(future, strikes, tau, sigmas) {
function (par) {
alpha <- par[1]
nu <- par[2]
rho <- par[3]
beta <- par[4]
sabr_sigs <- sapply(strikes, function(strike) {
sabrimpvol(fut = future, strike = strike, bizdays = tau,
alpha = alpha, beta = beta, rho = rho, nu = nu)
})
sum((sigmas - sabr_sigs)^2)
}
}
#' SABR model std. Deltas
#'
#' Finds and returns a dataframe with the strikes and implied
#' volatilites for the provided set of deltas or for the
#' default deltas.
#'
#' @param alpha a SABR model fit alpha parameter
#' @param beta a SABR model fit beta parameter
#' @param rho a SABR model fit rho parameter
#' @param nu a SABR model fit nu parameter
#' @param S underlying spot price at a given date
#' @param tau time in years from a given date to a maturity date
#' @param r interest free rate
#' @param q yield rate
#' @param future future price. Can be given instead of spot/rate/yield.
#' @param deltas deltas for which strikes and vols should be calculated. Defaults to usual
#' delta values.
#'
#' @return
#' Dataframe with columns delta, strike and impvol
#'
#' @export
sabr_deltas <- function(alpha, beta, rho, nu, S = NULL, tau, r = NULL,
future = NULL, q = NULL, deltas = DELTAS) {
eps <- 1e-8
if (is.null(future)) {
ft <- S * exp((r - q) * tau)
} else {
ft <- future
}
if (is.null(ft)) stop("Future price or spot, rate and yield must be provided.")
strk_sabr <- sapply(deltas/100, function(delta) {
uniroot(f=sabr_strikes, interval=c(1, ft*10), tol=eps,
t=tau, delta=delta, ft=ft, alpha=alpha,
beta=beta, rho=rho, nu=nu)$root
})
## multiroot não funciona bem, varx <- sabrimpvol(...) não usa vetor de deltas
# strk_sabr <- multiroot(func_ = sabr_strikes, interval = c(1, S*10), tolerance = eps,
# t = tau, delta = deltas, ft = ft, r = r, spot = S, alpha = alpha,
# beta = beta, rho = rho, nu = nu)$root
impvol_sabr <- sabrimpvol(fut = ft, strike = strk_sabr, bizdays = tau,
alpha = alpha, beta = beta, rho = rho, nu = nu)
data.frame(delta=deltas, strike=strk_sabr, impvol=impvol_sabr)
}
sabr_strikes <- function(x, delta, ft, t, alpha, beta, rho, nu) {
varx <- sabrimpvol(fut = ft, strike = x, bizdays = t,
alpha = alpha, beta = beta, rho = rho, nu = nu)
x - blackstrike(ft, delta, varx, t)
}
#' SVI Variance
#'
#' Calculates the variance as given by the SVI model for
#' the provided parameters.
#'
#' @param a 'a' in the SVI model (~~ instantaneous variance)
#' @param b 'b' in the SVI model (~~ speed of variance mean reversion)
#' @param m 'm' in the SVI model (~~ long term expected variance)
#' @param rho 'rho' in the SVI model (~~ correlation of variance and price)
#' @param x 'x' in the SVI model, given by ln(K/F) where K is an options strike price and F its underlying's future
#' @param sigma 'sigma' in the SVI model (~~ volatility of variance)
#'
#' @section Recycle rule:
#'
#' These arguments handle the recycle rule so vectors can be provided
#' and once those vectors differs in length the recycle rule is applied.
#'
#' @return
#' Volatility value as given by the SVI model
#'
#' @seealso
#' \code{\link{SVI_fit}} for SVI model fitting.
#' \code{\link{SVI_deltas}} for calculating strikes and imp. vols. for given deltas
#'
#' @export
SVI_var <- function(a, b, m, rho, x, sigma) {
a + b * (rho * (x - m) + sqrt((x - m) ^ 2 + sigma ^ 2))
}
#' SVI fit
#'
#' Fits a SVI model for the given smile as defined by strikes, futures
#' and their respective variances. Time to maturity should also be provided
#' for the non-arbitrage restriction. Returns the fitted parameters.
#'
#' @param strikes options strikes
#' @param future options underlying future If NULL will try to use spot and rate.
#' @param variances variances for the given strikes
#' @param ttm time to maturity in years
#' @param spot underlying spot priced, used in conjunction with rate if future price is not provided
#' @param rate risk free rate, used for calculating a future price in conjunction with spot if not provided
#'
#' @return
#' A list with the SVI model parameters a, b, x, m, rho and sigma.
#'
#' @seealso
#' \code{\link{SVI_var}} for SVI model variance calculation
#' \code{\link{SVI_deltas}} for calculating strikes and imp. vols. for given deltas
#'
#' @export
SVI_fit <- function(strikes, future = NULL, variances, ttm, spot = NULL, rate = NULL, initial_guess = NULL) {
if (is.null(future) && (is.null(spot) || is.null(rate))) stop("Either future price or spot and rate must be provided.")
if (is.null(future)) future <- spot * exp(rate * ttm)
x <- log(strikes/future)
tol <- length(variances) * min(variances) ^ 2
# SVI restrictions
# 0 <= a <= max(variances)
# 0 <= b <= 4/(ttm * (1 + abs(rho)))
# -1 <= rho <= 1
# min(x) <= m <= max(x)
# 0 < sigma < 10
# a, b, m, rho, sigma
eval_g0 <- function(par) {
c(par[2] - 4/(ttm * (1 + abs(par[4]))))
}
lo <- c(0, 0, min(x), -1, 1e-12)
hi <- c(max(c(variances, SVI.A)), 4/(ttm * (1 + 1e-4)), max(x), 1, 10)
if (is.null(initial_guess)) {
initial_guess <- c(a=SVI.A, b=SVI.B,
m=min(x), rho=SVI.RHO, sigma=SVI.SIGMA)
fit <- nloptr::nloptr(x0 = initial_guess,
eval_f = SVI_minsq(variances, x),
lb = lo,
ub = hi,
eval_g_ineq = eval_g0,
opts = list(algorithm = "NLOPT_GN_ISRES",
print_level = 0,
xtol_rel = 1e-4,
maxeval = 2000)
)
initial_guess <- fit$solution
}
fit <- nloptr::nloptr(x0 = initial_guess,
eval_f = SVI_minsq(variances, x),
lb = lo,
ub = hi,
eval_g_ineq = eval_g0,
opts = list(algorithm = "NLOPT_LN_COBYLA",
print_level = 0,
xtol_rel = 1e-4,
maxeval = 2000)
)
curr_param <- fit$solution
list(a = curr_param[1], b = curr_param[2], m = curr_param[3],
rho = curr_param[4], sigma = curr_param[5])
}
SVI_minsq <- function(variances, x) {
function (...) {
par <- c(...)
a <- par[1]
b <- par[2]
m <- par[3]
rho <- par[4]
sigma <- par[5]
SVI <- SVI_var(a, b, m, rho, x, sigma)
sum((variances - SVI)^2)
}
}
#' SVI fit std. Deltas
#'
#' Finds and returns a dataframe with the strikes and implied
#' volatilites for the provided set of deltas or for the
#' default deltas.
#'
#' @param a a SVI model fit a parameter
#' @param b a SVI model fit b parameter
#' @param m a SVI model fit m parameter
#' @param rho a SVI model fit rho parameter
#' @param sigma a SVI model fit sigma parameter
#' @param future future price
#' @param t time to maturity
#' @param deltas deltas for which strikes and vols should be calculated. Defaults to usual
#' delta values.
#'
#' @return
#' Dataframe with columns delta, strike and impvol
#'
#' @export
SVI_deltas <- function(a, b, m, rho, sigma, future, t, deltas = DELTAS) {
eps <- 1e-8
strk_svi <- sapply(deltas/100, function(delta) {
uniroot(f=svi_strikes, interval=c(1, future*10), tol=eps,
t = t, delta = delta, a = a, b = b,
m = m, rho = rho, sigma = sigma, ft = future)$root
})
impvol_svi <- sqrt(SVI_var(a, b, m, rho, log(strk_svi/future), sigma))
data.frame(delta=deltas, strike=strk_svi, impvol=impvol_svi)
}
svi_strikes <- function(x, delta, a, b, m, rho, sigma, ft, t) {
varx <- sqrt(SVI_var(a, b, m, rho, log(x/ft), sigma))
x - blackstrike(ft, delta, varx, t)
}
# calculate historical volatility for the period in
# ohlc_data using log returns
calc_hist_vol <- function (ohlc_data, year_biz_days) {
# ln(s_i+1 / s_i)
close_vals <- ohlc_data[, 1]
x <- sapply(2:nrow(ohlc_data), function(idx) {
log((close_vals[[idx]])/(close_vals[[idx - 1]]))
})
x_mean <- mean(x)
x_sum <- sum((x - x_mean)^2)
sqrt(x_sum * year_biz_days/(nrow(ohlc_data) - 2)) # no daily return for the first value
}
# Build volatility table for a given ohlc_data.
# Periods will be separated according to the number of
# months available in the data.
# Each row of the volatility table indicates when a period
# starts, and each column indicates the period length for
# which historical volatilities will be calculated.
build_vol_table <- function(ohlc_data) {
start_date <- as.Date(min(zoo::index(ohlc_data)))
end_date <- as.Date(max(zoo::index(ohlc_data)))
sample_size <- length(zoo::index(ohlc_data))
suppressWarnings(nmonth <- abs(interval(start_date, end_date)) %/% months(1))
days_by_period <- sample_size%/%nmonth
ymd_periods <- zoo::as.Date(sapply(0:(nmonth - 1), function(m) {zoo::index(ohlc_data[1 + m*days_by_period])}))
year_days <- bizdays(from=paste(year(start_date), '01', '01', sep='-'),
to=paste(year(start_date), '12', '31', sep='-'),
cal='Brazil/ANBIMA')
if(length(ymd_periods) > nmonth) ymd_periods <- ymd_periods[-length(ymd_periods)]
vol_table <- data.frame(row.names=ymd_periods)
for (m in 1:nmonth) {
vol_table[, as.character(m*days_by_period)] <- sapply(ymd_periods, function (st_date) {
st_date_index <- match(as.Date(st_date), zoo::index(ohlc_data))
ed_date_index <- st_date_index + days_by_period*m
if(ed_date_index <= (nmonth)*days_by_period + 1 && ed_date_index - st_date_index > 2) {
if (ed_date_index > sample_size) ed_date_index <- sample_size
tryCatch( {calc_hist_vol(ohlc_data[st_date_index:ed_date_index], year_days)}, error=function(e) {browser()})
} else NA
})
}
vol_table
}
# Calculate a confidence interval for a given alpha
confidence_interval <- function(alpha, n) {
lower <- sqrt((n - 1)/qchisq(p = alpha/2, df = n - 1))
higher <- sqrt((n - 1)/qchisq(p = (1 - alpha/2), df = n - 1))
c(lower, higher)
}
# Calculate implied volatilities for all options traded on ref_d
# Risk free rate is given by PRE curve at ref_d
# Implied volatility is calculated using bsm
calc_impl_vol <- function(options_data, ref_d) {
dividendYield <- 0
options <- options_data[!is.na(options_data$close), ]
options <- options[options$close != 0, ]
options <- options[options$ref_date==ref_d, ]
options <- transform(options, bznsdtexp = bizdays(ref_date, maturity_date, cal='Brazil/ANBIMA'))
options <- transform(options, dtexp = as.numeric(as.Date(maturity_date) - as.Date(ref_date)))
rate_curve <- get_curve('PRE', ref_d)
if(nrow(options) >= 1) {
impl_vol <- apply(options, 1, function(option) {
close <- as.numeric(option[['close']])
type <- option[['type']]
spot <- as.numeric(option[['spot']])
strike <- as.numeric(option[['strike']])
dtexp <- as.numeric(option[['dtexp']])
bzdtexp <- as.numeric(option[['bznsdtexp']])
riskFreeRate <- rate_curve[rate_curve$terms==dtexp, ]$value
i <- dtexp - 5
while (length(riskFreeRate) == 0) {
riskFreeRate <- rate_curve[rate_curve$terms==i, ]$value
i <- i + 1
}
ivol <- MyQuantTK::bsmimpvol(close, type, spot, strike, dtexp/365, riskFreeRate/100, dividendYield)
data.frame(bz_days_to_exp=bzdtexp, days_to_exp=dtexp, value=ivol, type=type)
})
impl_vol <- do.call(rbind, impl_vol)
} else {
impl_vol <- NULL
}
impl_vol
}
#' Volatility cone graph generator
#'
#' Generates a graph showing a volatility cone for the stock given
#' by ticker, calculated with data starting from historical_sample_start_date
#' to the most recent available date. It also calculates implied volatilities
#' for that stock's options that were traded on ref_date and plots it on
#' the same graph.
#' The volatility estimator for the cone is annualized standard deviation
#' of log returns. Its confidence interval is for a 5% alpha.
#' The implied volatilities are calculated using the bissection method
#' and the Black-Scholes model.
#'
#' @param ticker Stock to be graphed
#' @param ref_date Reference date for the options and impl. vol calculation
#' @param historical_sample_start_date Start date for the historical data to be used
#' to calculate the volatility cone
#'
#' @return
#' ggplot object containing the graph with the estimated volatility cone, its
#' confidence intervals and implied volatilities for options traded on ref_date
#'
#'
#' @export
gen_vol_cone <- function(ticker, ref_date, historical_sample_start_date, historical_sample_end_date=today()) {
options <- get_options(ticker, src='stock', from=ref_date, to=ref_date)
vals <- get_stocks(ticker, src='bloomberg', from=historical_sample_start_date, to=historical_sample_end_date)
vt <- build_vol_table(vals)
impl_vols <- calc_impl_vol(options, ref_date)
if(!is.null(impl_vols)) {
impl_vols$value <- impl_vols$value*100 # adjust to percentage values
# remove entries for which implied volatilities couldn't be found
# with the bissection method
impl_vols <- impl_vols[impl_vols$value!=0.000001, ]
}
# Adjust historical volatility data and create confidence interval
# with a 5% alpha
max_values <- apply(vt[, names(vt)], 2, max, na.rm = TRUE)*100
mean_values <- apply(vt[, names(vt)], 2, mean, na.rm = TRUE)*100
min_values <- apply(vt[, names(vt)], 2, min, na.rm = TRUE)*100
ci <- confidence_interval(0.05, length(vals))
maxcilow <- c(max_values*ci[1])
maxcihi <- c(max_values*ci[2])
meancilow <- c(mean_values*ci[1])
meancihi <- c(mean_values*ci[2])
mincilow <- c(min_values*ci[1])
mincihi <- c(min_values*ci[2])
cone <- melt(data.frame(days_to_exp=as.numeric(names(vt)),
max_values = max_values,
maxcilow = maxcilow,
maxcihi = maxcihi,
mean_values = mean_values,
meancilow = meancilow,
meancihi = meancihi,
min_values = min_values,
mincilow = mincilow,
mincihi = mincihi)
,id.vars = "days_to_exp")
len <- length(max_values)
cone$cgroup <- c(rep('Maximum', times=3*len), rep('Mean', times=3*len), rep('Minimum', times=3*len))
cone$type <- rep(c(rep('Estimate', times=len), rep('Confidence Interval', times=2*len)), times=3)
cone$type_detail <- rep(c(rep('Estimate', times=len), rep('Confidence Interval lo', times=len), rep('Confidence Interval hi', times=len)), times=3)
title <- paste(ticker, 'estimated vol. cone with impl. vols for', ref_date)
subtitle <- paste('Historical volatility calculated between', historical_sample_start_date, 'and', historical_sample_end_date)
plot <- ggplot(cone, aes(x=days_to_exp, y=value,
group=interaction(cgroup, type_detail),
colour=cgroup, linetype=type)) +
geom_line() +
geom_point(data = subset(cone, type == "Estimate")) +
scale_linetype_manual(values=c("dotted", "solid")) +
scale_x_continuous("Days to Expiry", breaks = c(0, as.numeric(names(vt)))) +
scale_y_continuous("Vol. (%)") +
coord_cartesian(ylim=c(max(c(min(mincihi - 5), 0)), max(maxcilow) + 5)) +
labs(shape='Option type', linetype="Line type", colour="Estimated historical vols.", title=title, subtitle=subtitle) +
theme(plot.title = element_text(hjust = 0.5))
if(!is.null(impl_vols)) {
plot <- plot + geom_point(data=impl_vols, aes(x=bz_days_to_exp, y=value, shape=type), inherit.aes=FALSE)
}
plot
}
#' volatility Class Coercing
#'
#' Generic function for coercing to \code{"volatility"}. There is a default
#' method and a method for \code{"variance"}.
#'
#' @param x an object whose class will determine the method to be dispatched.
#' Value is the volatility or variance.
#' @param ... further arguments to be passed to the next method, such as:
#' @param annualized for default method, TRUE if volatility value in x is annualized,
#' FALSE if volatility is daily.
#' @param dib days in one year, day count convention.
#'
#' @return
#' All values of x as volatility class, according to given daily or annualized type.
#'
#' @seealso
#' \code{\link{as.variance}}.
#'
#' @examples
#' as.volatility(1)
#' as.volatility(as.variance(.5))
#' as.volatility(0.2, TRUE, 252)
#' as.volatility(as.variance(10:30/100, TRUE))
#'
#' @name as.volatility
#' @export
as.volatility <- function(x, ...) UseMethod('as.volatility')
#' @rdname as.volatility
#' @export
as.volatility.default <- function(x, annualized=FALSE, dib=NULL, ...)
structure(x, annualized=annualized, dib=dib, class='volatility')
#' @rdname as.volatility
#' @export
as.volatility.variance <- function(x, ...) {
y <- unclass(x)
attributes(y) <- NULL
p <- if (attr(x, 'annualized')) 100 else 1
structure(sqrt(y) * p, annualized=attr(x, 'annualized'), dib=attr(x, 'dib'),
class='volatility')
}
#' variance Class Coercing
#'
#' Generic function for coercing to \code{"variance"}. There is a default
#' method and a method for \code{"volatility"}.
#'
#' @param x an object whose class will determine the method to be dispatched.
#' Value is the volatility or variance.
#' @param ... further arguments to be passed to the next method, such as:
#' @param annualized for default method, TRUE if variance value in x is annualized,
#' FALSE if variance is daily.
#' @param dib days in one year, day count convention.
#'
#' @return
#' All values of x as \code{"variance"} class, according to given daily or
#' annualized type.
#'
#' @seealso
#' \code{\link{as.volatility}}.
#'
#' @examples
#' as.variance(1)
#' as.variance(as.volatility(.5))
#' as.variance(0.2, TRUE, 252)
#' as.variance(as.volatility(100:300/100, TRUE))
#'
#' @name as.variance
#' @export
as.variance <- function(x, ...) UseMethod('as.variance')
#' @rdname as.variance
#' @export
as.variance.default <- function(x, annualized=FALSE, dib=NULL, ...)
structure(x, annualized=annualized, dib=dib, class='variance')
#' @rdname as.variance
#' @export
as.variance.volatility <- function(x, ...) {
y <- unclass(x)
attributes(y) <- NULL
p <- if (attr(x, 'annualized')) 100 else 1
structure((y / p) ^ 2, annualized=attr(x, 'annualized'), dib=attr(x, 'dib'),
class='variance')
}
#' Special print method for volatility class.
#'
#' Makes objects of class \code{"volatility"} be printed differently.
#'
#' @param x \code{"volatility"} class object to be printed.
#' @param ... further arguments to be passed to the next method.
#'
#' @return
#' invisible(x).
#'
#' @seealso
#' \code{\link{print.variance}}.
#'
#' @examples
#' a <- as.volatility(1)
#' print(a)
#' a <- annualize(a)
#' print(a)
#'
#' @export
print.volatility <- function(x, ...) {
cat(if (attr(x, "annualized")) "Annual" else "Daily",
"Volatility",
if (attr(x, "annualized")) paste("(%)", attr(x, "dib"), "days"),
"\n")
y <- unclass(x)
attributes(y) <- NULL
print(y)
invisible(x)
}
#' Special print method for variance class.
#'
#' Makes objects of class \code{"variance"} be printed differently.
#'
#' @param x \code{"variance"} class object to be printed.
#' @param ... further arguments to be passed to the next method.
#'
#' @return
#' invisible(x).
#'
#' @seealso
#' \code{\link{print.volatility}}.
#'
#' @examples
#' a <- as.variance(1)
#' print(a)
#' a <- annualize(a)
#' print(a)
#'
#' @export
print.variance <- function(x, ...) {
cat(if (attr(x, "annualized")) "Annual" else "Daily",
"Variance",
if (attr(x, "annualized")) paste(attr(x, "dib"), "days"),
"\n")
y <- unclass(x)
attributes(y) <- NULL
print(y)
invisible(x)
}
#' Annualize Volatility or Variance
#'
#' Transforms variance or volatility to annual. There is no default
#' method, so only classes \code{"volatility"}, \code{"variance"}, are supported.
#'
#' @param x an object whose class will determine the method to be dispatched.
#' Volatility or variance, will become annual.
#' @param ... further arguments to be passed to the next method.
#' @param dib days in one year, day count convention.
#'
#' @return
#' Same class with annual volatility or variance.
#'
#' @seealso
#' \code{\link{daily}} to transform variance or volatility to daily basis.
#'
#' @examples
#' annualize(as.variance(.5))
#' annualize(as.volatility(.25))
#'
#' @name annualize
#' @export
annualize <- function(x, ...) UseMethod('annualize')
#' @rdname annualize
#' @export
annualize.volatility <- function(x, dib=252, ...) {
y <- unclass(x)
attributes(y) <- NULL
if (attr(x, 'annualized'))
return(x)
else
y <- y * sqrt(dib) * 100
as.volatility(y, annualized=TRUE, dib=dib)
}
#' @rdname annualize
#' @export
annualize.variance <- function(x, dib=252, ...) {
y <- unclass(x)
attributes(y) <- NULL
if (attr(x, 'annualized'))
return(x)
else
y <- y * dib
as.variance(y, annualized=TRUE, dib=dib)
}
#' Daily Volatility or Variance
#'
#' Transforms variance or volatility to daily. There is no default
#' method, so only classes \code{"volatility"}, \code{"variance"}, are supported.
#'
#' @param x an object whose class will determine the method to be dispatched.
#' Volatility or variance, will become daily.
#' @param ... further arguments to be passed to the next method.
#'
#' @return
#' Same class with daily volatility or variance.
#'
#' @seealso
#' \code{\link{annualize}} to transform variance or volatility to annual basis.
#'
#' @examples
#' daily(annualize(as.variance(.5)))
#' daily(annualize(as.volatility(.25)))
#'
#' @name daily
#' @export
daily <- function(x, ...) UseMethod('daily')
#' @rdname daily
#' @export
daily.volatility <- function(x, ...) {
y <- unclass(x)
attributes(y) <- NULL
if (!attr(x, 'annualized'))
return(x)
else
y <- y / (sqrt(attr(x, 'dib')) * 100)
as.volatility(y, annualized=FALSE, dib=NULL)
}
#' @rdname daily
#' @export
daily.variance <- function(x, ...) {
y <- unclass(x)
attributes(y) <- NULL
if (! attr(x, 'annualized'))
return(x)
else
y <- y / attr(x, 'dib')
as.variance(y, annualized=FALSE, dib=NULL)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment