https://github.com/cran/fOptions
Tip revision: 9c58bb98a62fa8226a1f41c6fb1f98485b101250 authored by Rmetrics Core Team on 08 August 1977, 00:00:00 UTC
version 270.74
version 270.74
Tip revision: 9c58bb9
BasicAmericanOptions.R
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General
# Public License along with this library; if not, write to the
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA 02111-1307 USA
# Copyrights (C)
# for this R-port:
# 1999 - 2004, Diethelm Wuertz, GPL
# Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# info@rmetrics.org
# www.rmetrics.org
# for the code accessed (or partly included) from other R-ports:
# see R's copyright and license files
# for the code accessed (or partly included) from contributed R-ports
# and other sources
# see Rmetrics's copyright file
################################################################################
# FUNCTION: DESCRIPTION:
# RollGeskeWhaleyOption Roll-Geske-Whaley Calls on Dividend Paying Stocks
# BAWAmericanApproxOption Barone-Adesi and Whaley Approximation
# BSAmericanApproxOption Bjerksund and Stensland Approximation
################################################################################
RollGeskeWhaleyOption =
function(S, X, time1, Time2, r, D, sigma, title = NULL, description = NULL)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the option price of an American call on a stock
# paying a single dividend with specified time to divident
# payout. The option valuation formula derived by Roll, Geske
# and Whaley is used.
# References:
# Haug E.G., The Complete Guide to Option Pricing Formulas
# FUNCTION:
# Settings:
big = 100000000
eps = 1.0e-5
t1 = time1
T2 = Time2
# Compute:
Sx = S - D * exp(-r * t1)
if(D <= X * (1 - exp(-r*(T2-t1)))) {
result = GBSOption("c", Sx, X, T2, r, b=r, sigma)@price
cat("\nWarning: Not optimal to exercise\n")
return(result) }
ci = GBSOption("c", S, X, T2-t1, r, b=r, sigma)@price
HighS = S
while ( ci-HighS-D+X > 0 && HighS < big ) {
HighS = HighS * 2
ci = GBSOption("c", HighS, X, T2-t1, r, b=r, sigma)@price }
if(HighS > big) {
result = GBSOption("c", Sx, X, T2, r, b=r, sigma)@price
stop()}
LowS = 0
I = HighS * 0.5
ci = GBSOption("c", I, X, T2-t1, r, b=r, sigma)@price
# Search algorithm to find the critical stock price I
while ( abs(ci-I-D+X) > eps && HighS - LowS > eps ) {
if(ci-I-D+X < 0 ) { HighS = I }
else { LowS = I }
I = (HighS + LowS) / 2
ci = GBSOption("c", I, X, T2-t1, r, b=r, sigma)@price }
a1 = (log(Sx/X) + (r+sigma^2/2)*T2) / (sigma*sqrt(T2))
a2 = a1 - sigma*sqrt(T2)
b1 = (log(Sx/I) + (r+sigma^2/2)*t1) / (sigma*sqrt(t1))
b2 = b1 - sigma*sqrt(t1)
result = Sx*CND(b1) + Sx*CBND(a1,-b1,-sqrt(t1/T2)) -
X*exp(-r*T2)*CBND(a2,-b2,-sqrt(t1/T2)) -
(X-D)*exp(-r*t1)*CND(b2)
# Parameters:
# S, X, time1, Time2, r, D, sigma
param = list()
param$S = S
param$X = X
param$time1 = time1
param$Time2 = Time2
param$r = r
param$D = D
param$sigma = sigma
# Add title and description:
if(is.null(title)) title = "Roll Geske Whaley Option"
if(is.null(description)) description = as.character(date())
# Return Value:
new("fOPTION",
call = match.call(),
parameters = param,
price = result,
title = title,
description = description
)
}
# ******************************************************************************
BAWAmericanApproxOption =
function(TypeFlag = c("c", "p"), S, X, Time, r, b, sigma, title = NULL,
description = NULL)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the option price of an American call or put
# option on an underlying asset for a given cost-of-carry rate.
# The quadratic approximation method by Barone-Adesi and
# Whaley is used.
# References:
# Haug E.G., The Complete Guide to Option Pricing Formulas
# FUNCTION:
# Settings:
TypeFlag = TypeFlag[1]
# Compute:
if(TypeFlag == "c") {
result = .BAWAmCallApproxOption(S, X, Time, r, b, sigma) }
if(TypeFlag == "p") {
result = .BAWAmPutApproxOption(S, X, Time, r, b, sigma) }
# Parameters:
# TypeFlag = c("c", "p"), S, X, Time, r, b, sigma
param = list()
param$TypeFlag = TypeFlag
param$S = S
param$X = X
param$Time = Time
param$r = r
param$b = b
param$sigma = sigma
# Add title and description:
if(is.null(title)) title = "BAW American Approximated Option"
if(is.null(description)) description = as.character(date())
# Return Value:
new("fOPTION",
call = match.call(),
parameters = param,
price = result,
title = title,
description = description
)
}
.BAWAmCallApproxOption <-
function(S, X, Time, r, b, sigma)
{
# Internal Function - The Call:
# Compute:
if(b >= r) {
result = GBSOption("c", S, X, Time, r, b, sigma)@price }
else {
Sk = .bawKc(X, Time, r, b, sigma)
n = 2*b/sigma^2
K = 2*r/(sigma^2*(1-exp(-r*Time)))
d1 = (log(Sk/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
Q2 = (-(n-1)+sqrt((n-1)^2+4*K))/2
a2 = (Sk/Q2)*(1-exp((b-r)*Time)*CND(d1))
if(S < Sk) {
result = GBSOption("c", S, X, Time, r, b, sigma)@price +
a2*(S/Sk)^Q2
} else {
result = S-X
}
}
# Return Value:
result
}
.bawKc <-
function(X, Time, r, b, sigma)
{
# Newton Raphson algorithm to solve for the critical commodity
# price for a Call.
# Calculation of seed value, Si
n = 2*b/sigma^2
m = 2*r/sigma^2
q2u = (-(n-1)+sqrt((n-1)^2+4*m))/2
Su = X/(1-1/q2u)
h2 = -(b*Time+2*sigma*sqrt(Time))*X/(Su-X)
Si = X+(Su-X)*(1-exp(h2))
K = 2*r/(sigma^2*(1-exp(-r*Time)))
d1 = (log(Si/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
Q2 = (-(n-1)+sqrt((n-1)^2+4*K))/2
LHS = Si-X
RHS = GBSOption("c", Si, X, Time, r, b, sigma)@price +
(1-exp((b-r)*Time)*CND(d1))*Si/Q2
bi = exp((b-r)*Time)*CND(d1)*(1-1/Q2) +
(1-exp((b-r)*Time)*CND(d1)/(sigma*sqrt(Time)))/Q2
E = 0.000001
# Newton Raphson algorithm for finding critical price Si
while (abs(LHS-RHS)/X > E) {
Si = (X+RHS-bi*Si)/(1-bi)
d1 = (log(Si/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
LHS = Si-X
RHS = GBSOption("c", Si, X, Time, r, b, sigma)@price +
(1-exp((b-r)*Time)*CND(d1))*Si/Q2
bi = exp((b-r)*Time)*CND(d1)*(1-1/Q2) +
( 1-exp((b-r)*Time)*CND(d1)/(sigma*sqrt(Time)))/Q2 }
# Return Value:
Si
}
.BAWAmPutApproxOption <-
function(S, X, Time, r, b, sigma)
{
# Internal Function - The Put:
# Compute:
Sk = .bawKp(X, Time, r, b, sigma)
n = 2*b/sigma^2
K = 2*r/(sigma^2*(1-exp(-r*Time)))
d1 = (log(Sk/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
Q1 = (-(n-1)-sqrt((n-1)^2+4*K))/2
a1 = -(Sk/Q1)*(1-exp((b-r)*Time)*CND(-d1))
if(S > Sk) {
result = GBSOption("p", S, X, Time, r, b, sigma)@price + a1*(S/Sk)^Q1
} else {
result = X-S
}
# Return Value:
result
}
.bawKp <-
function(X, Time, r, b, sigma)
{
# Internal Function - used for the Put:
# Newton Raphson algorithm to solve for the critical commodity
# price for a Put.
# Calculation of seed value, Si
n = 2*b/sigma^2
m = 2*r/sigma^2
q1u = (-(n-1)-sqrt((n-1)^2+4*m))/2
Su = X/(1-1/q1u)
h1 = (b*Time-2*sigma*sqrt(Time))*X/(X-Su)
Si = Su+(X-Su)*exp(h1)
K = 2*r/(sigma^2*(1-exp(-r*Time)))
d1 = (log(Si/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
Q1 = (-(n-1)-sqrt((n-1)^2+4*K))/2
LHS = X-Si
RHS = GBSOption("p", Si, X, Time, r, b, sigma)@price -
(1-exp((b-r)*Time)*CND(-d1))*Si/Q1
bi = -exp((b-r)*Time)*CND(-d1)*(1-1/Q1) -
(1+exp((b-r)*Time)*CND(-d1)/(sigma*sqrt(Time)))/Q1
E = 0.000001
# Newton Raphson algorithm for finding critical price Si
while (abs(LHS-RHS)/X > E ) {
Si = (X-RHS+bi*Si)/(1+bi)
d1 = (log(Si/X)+(b+sigma^2/2)*Time)/(sigma*sqrt(Time))
LHS = X-Si
RHS = GBSOption("p", Si, X, Time, r, b, sigma)@price -
(1-exp((b-r)*Time)*CND(-d1))*Si/Q1
bi = -exp((b-r)*Time)*CND(-d1)*(1-1/Q1) -
(1+exp((b-r)*Time)*CND(-d1)/(sigma*sqrt(Time)))/Q1 }
# Return Value:
Si
}
# ------------------------------------------------------------------------------
BSAmericanApproxOption =
function(TypeFlag = c("c", "p"), S, X, Time, r, b, sigma, title = NULL,
description = NULL)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the option price of an American call or
# put option stocks, futures, and currencies. The
# approximation method by Bjerksund and Stensland is used.
# References:
# Haug E.G., The Complete Guide to Option Pricing Formulas
# FUNCTION:
# Settings:
TypeFlag = TypeFlag[1]
# The Bjerksund and Stensland (1993) American approximation:
if(TypeFlag == "c") {
result = .BSAmericanCallApprox(S, X, Time, r, b, sigma) }
if(TypeFlag == "p") {
# Use the Bjerksund and Stensland put-call transformation
result = .BSAmericanCallApprox(X, S, Time, r - b, -b, sigma) }
# Parameters:
# TypeFlag = c("c", "p"), S, X, Time, r, b, sigma
param = list()
param$TypeFlag = TypeFlag
param$S = S
param$X = X
param$Time = Time
param$r = r
param$b = b
param$sigma = sigma
if(!is.na(result$TriggerPrice)) param$TrigerPrice = result$TriggerPrice
# Add title and description:
if(is.null(title)) title = "BS American Approximated Option"
if(is.null(description)) description = as.character(date())
# Return Value:
new("fOPTION",
call = match.call(),
parameters = param,
price = result$Premium,
title = title,
description = description
)
}
.BSAmericanCallApprox <-
function(S, X, Time, r, b, sigma)
{
# Call Approximation:
if(b >= r) {
# Never optimal to exersice before maturity
result = list(
Premium = GBSOption("c", S, X, Time, r, b, sigma)@price,
TriggerPrice = NA)
} else {
Beta = (1/2 - b/sigma^2) + sqrt((b/sigma^2 - 1/2)^2 + 2*r/sigma^2)
BInfinity = Beta/(Beta-1) * X
B0 = max(X, r/(r-b) * X)
ht = -(b*Time + 2*sigma*sqrt(Time)) * B0/(BInfinity-B0)
# Trigger Price I:
I = B0 + (BInfinity-B0) * (1 - exp(ht))
alpha = (I-X) * I^(-Beta)
if(S >= I) {
result = list(
Premium = S-X,
TriggerPrice = I) }
else {
result = list(
Premium = alpha*S^Beta - alpha*.bsPhi(S,Time,Beta,I,I,r,b,sigma) +
.bsPhi(S,Time,1,I,I,r,b,sigma) - .bsPhi(S,Time,1,X,I,r,b,sigma) -
X*.bsPhi(S,Time,0,I,I,r,b,sigma) + X*.bsPhi(S,Time,0,X,I,r,b,sigma),
TriggerPrice = I) } }
result}
.bsPhi <-
function(S, Time, gamma, H, I, r, b, sigma)
{
# Utility function phi:
lambda = (-r + gamma*b + 0.5*gamma * (gamma-1)*sigma^2) * Time
d = -(log(S/H) + (b + (gamma-0.5)*sigma^2)*Time) /
(sigma*sqrt(Time))
kappa = 2 * b / (sigma^2) + (2*gamma - 1)
result = exp(lambda)*S^gamma *
(CND(d)-(I/S)^kappa*CND(d-2*log(I/S)/(sigma*sqrt(Time))))
# Return Value:
result
}
################################################################################