# 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 Description. 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 # 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 ################################################################################ # MOMENT MATCHING: DESCRIPTION: # MomentMatchedAsianOption Valuate moment matched option prices # .LevyTurnbullWakemanAsianOption Log-Normal Approximation # .MilevskyPosnerAsianOption Reciprocal-Gamma Approximation # .PosnerMilevskyAsianOption Johnson Type I Approximation # MomentMatchedAsianDensity Valuate moment matched option densities # .LevyTurnbullWakemanAsianDensity Log-Normal Approximation # .MilevskyPosnerAsianDensity Reciprocal-Gamma Approximation # .PosnerMilevskyAsianDensity Johnson Type I Approximation # GRAM CHARLIER SERIES EXPANSION: DESCRIPTION: # GramCharlierAsianOption Calculate Gram-Charlier option prices # .GramCharlierAsianDensity NA # STATE SPACE MOMENTS: DESCRIPTION: # AsianOptionMoments Methods to calculate Asian Moments # .DufresneAsianOptionMoments Moments from Dufresne's Formula # .AbrahamsonAsianOptionMoments Moments from Abrahamson's Formula # .TurnbullWakemanAsianOptionMoments First 2 Moments from Turnbull-Wakeman # .TolmatzAsianOptionMoments Asymptotic Behavior after Tolmatz # STATE SPACE DENSITIES: DESCRIPTION: # StateSpaceAsianDensity NA # .Schroeder1AsianDensity NA # .Schroeder2AsianDensity NA # .Yor1AsianDensity NA # .Yor2AsianDensity NA # .TolmatzAsianDensity NA # .TolmatzAsianProbability NA # PARTIAL DIFFERENTIAL EQUATIONS: DESCRIPTION: # PDEAsianOption PDE Asian Option Pricing # .ZhangAsianOption Asian option price by Zhang's 1D PDE # ZhangApproximateAsianOption # .VecerAsianOption Asian option price by Vecer's 1D PDE # LAPLACE INVERSION: DESCRIPTION: # GemanYorAsianOption Asian option price by Laplace Inversion # gGemanYor Function to be Laplace inverted # SPECTRAL EXPANSION: DESCRIPTION: # LinetzkyAsianOption Asian option price by Spectral Expansion # gLinetzky Function to be integrated # BOUNDS ON OPTION PRICES: DESCRIPTION: # BoundsOnAsianOption Lower and upper bonds on Asian calls # CurranThompsonAsianOption From Thompson's continuous limit # RogerShiThompsonAsianOption From Thompson's single integral formula # ThompsonAsianOption Thompson's upper bound # SYMMETRY RELATIONS: DESCRIPTION: # CallPutParityAsianOption Call-Put parity Relation # WithDividendsAsianOption Adds dividends to Asian Option Formula # TABULATED RESULTS: DESCRIPTION: # FuMadanWangTable Table from Fu, Madan and Wang's paper # FusaiTaglianiTable Table from Fusai und tagliani's paper # GemanTable Table from Geman's paper # LinetzkyTable Table from Linetzky's paper # ZhangTable Table from Zhang's paper # ZhangLongTable Long Table from Zhang's paper # ZhangShortTable Short Table from Zhang's paper ################################################################################ ################################################################################ # MOMENT MATCHING: MomentMatchedAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30, table = NA, method = c("LN", "RG", "JI")) { # A function implemented by Diethelm Wuertz # Description: # Calculates price for Asian options based on moment matching # LN: Levy-Turnbull-Wakeman Log-Normal Approximation # RG: Milevsky-Posner Reciprocal-Gamma Approximation # JI: Posner-Milevski Johnson Type I Approximation # FUNCTION: # Set Default TypeFlag and Method, if no other is selected: TypeFlag = TypeFlag[1] method = method[1] # Test for Table: if (is.data.frame(table)) { S = table[,1] X = table[,2] Time = table[,3] r = table[,4] sigma = table[,5] } call = rep(0, length=length(S)) # Log-Normal Approximation: if (method == "LN") { for ( i in 1:length(S) ) { moments = masian(Time = Time[i], r = r[i], sigma = sigma[i])$rawMoments moments = moments / Time[i]^(1:4) meanlog = ( 2*log(moments[1]) - log(moments[2])/2 ) sdlog = ( sqrt ( log(moments[2]) - 2*log(moments[1]) ) ) d2 = ( -log(X[i]/S[i]) + meanlog*Time[i] ) / ( sdlog*sqrt(Time[i]) ) d1 = d2 + sdlog*sqrt(Time[i]) call[i] = moments[1]*pnorm(d1)-(X[i]/S[i])*pnorm(d2) } } # Reciprocal-Gamma Approximation: if (method == "RG") { for ( i in 1:length(S) ) { moments = masian(Time = Time[i], r = r[i], sigma = sigma[i])$rawMoments moments = moments / Time[i]^(1:4) alpha = (2*moments[2] - moments[1]^2) / (moments[2] - moments[1]^2) beta = (moments[2] - moments[1]^2) / (moments[1]*moments[2]) call[i] = moments[1]*pgamma(S[i]/X[i], shape=alpha-1, scale=beta) - (X[i]/S[i])*pgamma(S[i]/X[i], shape=alpha, scale=beta) } } # Johnson-Type-I Approximation: if (method == "JI") { for ( i in 1:length(S) ) { cmoments = masian(Time = Time[i], r = r[i], sigma = sigma[i])$centralMoments cmoments = cmoments / Time[i]^(1:4) mu = cmoments[1] varsigma = sqrt(cmoments[2]) eta = cmoments[3] / varsigma^3 omega = 0.5 * ( 8 + 4*eta^2 + 4*sqrt(4*eta^2+eta^4) )^(1/3) omega = omega + 1/omega - 1 b = 1 / sqrt(log(omega)) a = 0.5 * b * log(omega*(omega-1)/varsigma^2) d = sign(eta) c = d*mu - exp( (0.5/b-a)/b ) Q = a + b*log((X[i]/S[i]-c)/d) call[i] = mu - X[i]/S[i] + (X[i]/S[i] - c) * pnorm( Q ) - d * exp ( (1-2*a*b)/(2*b^2) ) * pnorm ( Q - 1/b ) } } # Call Price: Price = S* exp(-r*Time) * call # Put Price: if (TypeFlag == "p") { Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X Price = Price - Parity } # Return Value: option = list( price = Price, call = match.call() ) class(option) = "option" option } # ------------------------------------------------------------------------------ .LevyTurnbullWakemanAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Return Value: MomentMatchedAsianOption(TypeFlag[1], S, X, Time, r, sigma, method = "LN") } # ------------------------------------------------------------------------------ .MilevskyPosnerAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Return Value: MomentMatchedAsianOption(TypeFlag[1], S, X, Time, r, sigma, method = "RG") } # ------------------------------------------------------------------------------ .PosnerMilevskyAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Return Value: MomentMatchedAsianOption(TypeFlag[1], S, X, Time, r, sigma, method = "JI") } # ------------------------------------------------------------------------------ MomentMatchedAsianDensity = function(x, Time = 1, r = 0.09, sigma = 0.30, method = c("LN", "RG", "JI")) { # A function implemented by Diethelm Wuertz # Description: # Calculates price for Asian options based on moment matching # LN: Levy-Turnbull-Wakeman Log-Normal Approximation # RG: Milevsky-Posner Reciprocal-Gamma Approximation # JI: Posner-Milevski Johnson Type I Approximation # FUNCTION: # Set Default Method, if no other is selected: method = method[1] # Log-Normal Approximation: if (method == "LN") { moments = masian(Time = Time, r = r, sigma = sigma)$rawMoments moments = moments / Time^(1:4) meanlog = 2*log(moments[1]) - log(moments[2])/2 sdlog = sqrt ( log(moments[2]) - 2*log(moments[1]) ) density = dlnorm(x = x, meanlog = meanlog, sdlog = sdlog) } # Reciprocal-Gamma Approximation: if (method == "RG") { moments = masian(Time = Time, r = r, sigma = sigma)$rawMoments moments = moments / Time^(1:4) alpha = (2*moments[2] - moments[1]^2) / (moments[2] - moments[1]^2) beta = (moments[2] - moments[1]^2) / (moments[1]*moments[2]) density = drgam(x = x, alpha = alpha, beta = beta) } # Johnson Type I Approximation: if (method == "JI") { cmoments = masian(Time = Time, r = r, sigma = sigma)$centralMoments cmoments = cmoments / Time^(1:4) mu = cmoments[1] varsigma = sqrt(cmoments[2]) eta = cmoments[3] / varsigma^3 omega = 0.5 * ( 8 + 4*eta^2 + 4*sqrt(4*eta^2+eta^4) )^(1/3) omega = omega + 1/omega - 1 b = 1 / sqrt(log(omega)) a = 0.5 * b * log(omega*(omega-1)/varsigma^2) d = sign(eta) c = d*mu - exp( (0.5/b-a)/b ) density = djohnson(x = x, a = a, b = b, c = c, d = d) } # Return Value: density } # ------------------------------------------------------------------------------ .LevyTurnbullWakemanAsianDensity = function(x, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Return Value: MomentMatchedAsianDensity(x, Time, r, sigma, method = "LN") } # ------------------------------------------------------------------------------ .MilevskyPosnerAsianDensity = function(x, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Return Value: MomentMatchedAsianDensity(x, Time, r, sigma, method = "RG") } # ------------------------------------------------------------------------------ .PosnerMilevskyAsianDensity = function(x, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Return Value: MomentMatchedAsianDensity(x, Time, r, sigma, method = "JI") } ################################################################################ # GRAM CHARLIER: GramCharlierAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30, table = NA, method = c("LN", "RG", "JI")) { # A function implemented by Diethelm Wuertz # Description: # Calculate arithmetic Asian options price using # Gram Charlier Statistical Series Expansion around # LN: Log-Normal Distribution # RG: Reciprocal-Gamma Distribution # JI: Johnson-Type-I Distribution # FUNCTION: # Select Method: TypeFlag = TypeFlag[1] method = method[1] # Test for Table: if (is.data.frame(table)) { S = table[,1] X = table[,2] Time = table[,3] r = table[,4] sigma = table[,5] } # Calculate Price: Price = MomentMatchedAsianOption("c", S = S, X = X, Time = Time, r = r, sigma = sigma, method=method)$price gc3 = gc4 = rep(0, length(Price)) # Log-Normal Approximation: if (method == "LN") { for ( i in 1:length(S) ) { moments = masian(Time[i], r[i], sigma[i])$rawMoments/Time[i]^(1:4) meanlog = 2*log(moments[1]) - log(moments[2])/2 sdlog = sqrt ( log(moments[2]) - 2*log(moments[1]) ) asian.cm = masian(Time[i], r[i], sigma[i])$centralMoments/Time[i]^(1:4) lnorm.cm = mlognorm(meanlog, sdlog)$centralMoments/Time[i]^(1:4) kappa = (asian.cm-lnorm.cm) gc3[i] = kappa[3]*dlognorm(X[i]/S[i], meanlog, sdlog, deriv=1)/6 gc4[i] = kappa[4]*dlognorm(X[i]/S[i], meanlog, sdlog, deriv=2)/24 } } # Reciprocal-Gamma Approximation: if (method == "RG" ) { for ( i in 1:length(S) ) { moments = masian(Time[i], r[i], sigma[i])$rawMoments/Time[i]^(1:4) alpha = (2*moments[2] - moments[1]^2) / (moments[2] - moments[1]^2) beta = (moments[2] - moments[1]^2) / (moments[1]*moments[2]) asian.cm = masian(Time[i], r[i], sigma[i])$centralMoments/Time[i]^(1:4) rgam.cm = mrgam(alpha, beta)$centralMoments/Time[i]^(1:4) kappa = (asian.cm-rgam.cm) gc3[i] = kappa[3]*drgam(X[i]/S[i], alpha, beta, deriv = 1)/6 gc4[i] = kappa[4]*drgam(X[i]/S[i], alpha, beta, deriv = 2)/24 } } # Johnson-Type-I Approximation: if (method == "JI" ) { for ( i in 1:length(S) ) { cmoments = masian(Time = Time[i], r = r[i], sigma = sigma[i])$centralMoments cmoments = cmoments / Time[i]^(1:4) mu = cmoments[1] varsigma = sqrt(cmoments[2]) eta = cmoments[3] / varsigma^3 omega = 0.5 * ( 8 + 4*eta^2 + 4*sqrt(4*eta^2+eta^4) )^(1/3) omega = omega + 1/omega - 1 b = 1 / sqrt(log(omega)) a = 0.5 * b * log(omega*(omega-1)/varsigma^2) d = sign(eta) c = d*mu - exp( (0.5/b-a)/b ) asian.cm = cmoments johnson.cm = cmoments skewness = sqrt((omega-1)*(omega+2)^2) kurtosis = omega^4 + 2*omega^3 + 3* omega^2 - 3 johnson.cm[3] = skewness * varsigma^3 johnson.cm[4] = kurtosis * varsigma^4 kappa = (asian.cm-johnson.cm) gc3[i] = kappa[3]*djohnson(X[i]/S[i], a, b, c, d, deriv = 1)/6 gc4[i] = kappa[4]*djohnson(X[i]/S[i], a, b, c, d, deriv = 2)/24 } } # Gram-Charlier Approximated Call Price: Price = Price - S * exp(-r*Time) * (gc3-gc4) # Put Price: if (TypeFlag == "p") { Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X Price = Price - Parity } # Return Value: option = list( price = Price, call = match.call() ) class(option) = "option" option } .GramCharlierAsianDensity = function(Time = 1, r = 0.09, sigma = 0.30, method = c("LN", "RG", "JI")) { # A function ported by Diethelm Wuertz # Return Value: NA } ################################################################################ # STATE SPACE MOMENTS: AsianOptionMoments = function(M = 4, Time = 1, r = 0.045, sigma = 0.30, log = FALSE, method = c("A", "D", "TW", "T")) { # A function implemented by Diethelm Wuertz # Description: # Calculates Asian Moments using several approaches: # A - Moments from Abrahamson's Formula # D - Moments from Dufresne's Formula # TW - First 2 Moments from Turnbull-Wakeman # T - Asymptotic Behavior after Tolmatz # FUNCTION: # Settings: method = method[1] result = NA # Abrahamson Formula: if (method == "A") result = .AbrahamsonAsianOptionMoments(M = M, Time = Time, r = r, sigma = sigma) # Dufresne Formula: if (method == "D") result = .DufresneAsianOptionMoments(M = M, Time = Time, r = r, sigma = sigma) # Tolmatz Formula - Asymptotic Behavior: if (method == "T") result = .TolmatzAsianOptionMoments(M = M, Time = Time, r = r, sigma = sigma, log = log) # Turnbull Wakeman - Explicit 1st and Second Moment: if (method == "TW") result = .TurnbullWakemanAsianOptionMoments(M = M, Time = Time, r = r, sigma = sigma) # Return Value: result } # ------------------------------------------------------------------------------ .DufresneAsianOptionMoments = function(M = 4, Time = 1, r = 0.045, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Description: # Calculates Moments of Asian Options Density # according to the formula of Dufresne-GemanYor # FUNCTION: # Calculates: E[(A_\tau^{(\nu)})^n] moments = function (M, tau, nu) { d = function(j, n, beta) { d = 2^n for (i in 0:n) if (i != j) d = d / ( (beta+j)^2 - (beta+i)^2 ) d } moments = rep(0, length=M) for (n in 1:M) { moments[n] = 0 for (j in 0:n) moments[n] = moments[n] + d(j, n, nu/2)*exp(2*(j^2+j*nu)*tau) moments[n] = prod(1:n) * moments[n] / (2^(2*n)) } moments } # Calculate for: tau = sigma^2*Time/4 nu = 2*r/sigma^2-1 # Return Value: (4/sigma^2)^(1:M) * moments(M, tau, nu) } # ------------------------------------------------------------------------------ .AbrahamsonAsianOptionMoments = function (M = 4, Time = 1, r = 0.045, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Description: # Calculates Moments of Asian Options Density # according to the formula of Abrahamson # FUNCTION: # Calculates: E[(A_\tau^{(\nu)})^n] moments = function (M, tau, nu) { moments = rep(0, times = M) for (N in 1:M) { d = c = 2 * ( (1:N)^2 + (1:N)*nu ) for (m in 1:N) { for (j in 1:N) { if (j!= m) d[m] = d[m]*(c[m]-c[j]) } d[m] = exp(c[m]*tau) / d[m] } moments[N] = prod(1:N) * ( sum(d) + (-1)^N/prod(c) ) } moments } # Calculate for: tau = sigma^2*Time/4 nu = 2*r/sigma^2-1 # Return Value: (4/sigma^2)^(1:M) * moments(M, tau, nu) } # ------------------------------------------------------------------------------ .TurnbullWakemanAsianOptionMoments = function (M = 2, Time = 1, r = 0.045, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Description: # Calculates the first two moments as derived explicitly # by Turnbull and Wakeman. It can serve as a test for # other implementations. # Note: # Maximum M is 2! # FUNCTION: # Moments: moments = rep(NA, times = M) if (M == 1 || M == 2) moments[1] = (exp(r*Time)-1)/(r*Time) if (M == 2) moments[2] = 2*exp((2*r+sigma^2)*Time)/ ((r+sigma^2)*(2*r+sigma^2)*Time^2) + (2/(r*Time^2)) * ( 1/(2*r+sigma^2) - exp(r*Time)/(r+sigma^2) ) # Return Value: moments } # ------------------------------------------------------------------------------ .TolmatzAsianOptionMoments = function (M = 100, Time = 1, r = 0.045, sigma = 0.30, log = FALSE) { # A function implemented by Diethelm Wuertz # Description: # Calculates Asymptotic Moments of Asian Options Density # according to the formula of Tolmatz for nu=0 and Wuertz # for nu different from zero - Log returns can be selected # FUNCTION: # Calculates: log { E[(A_\tau^{(\nu)})^n] } moments = function (M, tau, nu=0) { moments = rep(0, times=M) M = 1:M log.moments = -M*log(2) + lgamma(nu+M) - lgamma(nu+2*M) + 2*(M^2+M*nu)*tau log.moments } # Calculate for: tau = sigma^2*Time/4 nu = 2*r/sigma^2-1 # Return Value: moments = (1:M)*log(4/sigma^2) + moments(M, tau, nu) # Return value: if (!log) moments = exp(moments) moments } ################################################################################ # ASIAN DENSITY: # STATE SPACE DENSITIES: DESCRIPTION: # StateSpaceAsianDensity # .Schroeder1AsianDensity S1 # .Schroeder2AsianDensity S2 # .Yor1AsianDensity Y1 # .Yor2AsianDensity Y2 # .TolmatzAsianDensity T # .TolmatzAsianProbability ################################################################################ # PDE SOLVER: ZhangAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30, table = NA, correction = TRUE, nint = 800, eps = 1.0e-8, dt = 1.0e-10) { # A function implemented by Diethelm Wuertz # Description: # Valuates Asian options by Solving Zhang's one # dimensional Partial Differential Equations # Source: # For the Fortran Routine: # TOMS ... # FUNCTION: # Settings: TypeFlag = TypeFlag[1] # Test for Table: if (is.data.frame(table)) { S = table[, 1] X = table[, 2] Time = table[, 3] r = table[, 4] sigma = table[, 5] } Price = rep(0, times = length(S)) # Set Model Identifier: modsel = 2 # Option Parameters: if (TypeFlag == "c") z = +1 if (TypeFlag == "p") z = -1 # PDE Parameters - Do not change: T0 = 0; Tout = 1 np = 0; Price.by.S = 0 mf = 12 npde = 1; kord = 4; ncc = 2; maxder = 5 # Fill Working Arrays: xbkpt = rep(0, times = nint+1) length.work = kord+npde*(4+9*npde)+(kord+(nint-1)*(kord-ncc)) * (3*kord+2+npde*(3*(kord-1)*npde+maxder+4)) work = rep(0, times = length.work) length.iwork = (npde+1)*(kord+(nint-1)*(kord-ncc)) iwork = rep(0, times = length.iwork) # Compute Prices: for ( i in 1:length(S) ) { result = .Fortran("asianval", as.double(z), as.double(S[i]), as.double(X[i]), as.double(X[i]), as.double(X[i]), as.double(Time[i]), as.double(r[i]), as.double(sigma[i]), as.double(T0), as.double(Tout), as.double(eps), as.double(dt), as.double(Price.by.S), as.integer(np), as.integer(modsel), as.integer(mf), as.integer(npde), as.integer(kord), as.integer(nint), as.integer(ncc), as.integer(maxder), as.double(X[i]/S[i]), as.double(xbkpt), as.double(work), as.integer(iwork), PACKAGE = "fAsianOptions" ) Price[i] = result[[13]]*S[i] } # ? Price = Price + ZhangApproximateAsianOption(TypeFlag, S, X, Time, r, sigma, table) # Return Value: Price } ZhangApproximateAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30, table = NA) { # Settings: TypeFlag = TypeFlag[1] # Test for Table: if (is.data.frame(table)) { S = table[, 1] X = table[, 2] Time = table[, 3] r = table[, 4] sigma = table[, 5] } # Compute: I = 0 xi = (Time*X-I)*exp(-r*Time)/S - (1-exp(-r*Time))/r eta = (sigma^2/(4*r^3)) * (-3 + 2*r*Time + 4*exp(-r*Time) - exp(-2*r*Time)) # Call: price = (S/Time) * ( -xi * pnorm(-xi/sqrt(2*eta)) + sqrt(eta/pi)*exp(-xi^2/(4*eta)) ) if (TypeFlag == "c") { ans = price } else { ans = CallPutParityAsianOption(TypeFlag = "c", price, S, X, Time, sigma, r, table = table) } # Return Value: ans } # ------------------------------------------------------------------------------ VecerAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30, table = NA, nint = 800, eps = 1.0e-8, dt = 1.0e-10) { # A function implemented by Diethelm Wuertz # Description: # Valuates Asian options by Solving Vecer's one # dimensional Partial Differential Equations # FUNCTION: # Vecer's PDE modeled by: modsel == 1 # SUBROUTINE ASIANVAL( # ZZ, SS, XS, XSMIN, XSMAX, TIME, RR, SIGMA, # T0, TOUT, EPS, DT, PRICEBYS, NP, MODSEL, # MF1, NPDE1, KORD1, MX1, NCC1, MAXDER1, # XBYS, XBKPT, WORK, IWORK) # Settings: TypeFlag = TypeFlag[1] # Test for Table: if (is.data.frame(table)) { S = table[, 1] X = table[, 2] Time = table[, 3] r = table[, 4] sigma = table[, 5] } Price = rep(0, times = length(S)) # Set Model Identifier: modsel = 1 # Option Parameters: if (TypeFlag == "c") z = +1 if (TypeFlag == "p") z = -1 # PDE Parameters - Do not change: T0 = 0 Tout = 1 np = 0 Price.by.S = 0 mf = 12 npde = 1 kord = 4 ncc = 2 maxder = 5 # Fill Working Arrays: xbkpt = rep(0, times = nint+1) length.work = kord+npde*(4+9*npde)+(kord+(nint-1)*(kord-ncc)) * (3*kord+2+npde*(3*(kord-1)*npde+maxder+4)) work = rep(0, times = length.work) length.iwork = (npde+1)*(kord+(nint-1)*(kord-ncc)) iwork = rep(0, times = length.iwork) # Compute Prices: for ( i in 1:length(S) ) { result = .Fortran("asianval", as.double(z), as.double(S[i]), as.double(X[i]), as.double(X[i]), as.double(X[i]), as.double(Time[i]), as.double(r[i]), as.double(sigma[i]), as.double(T0), as.double(Tout), as.double(eps), as.double(dt), as.double(Price.by.S), as.integer(np), as.integer(modsel), as.integer(mf), as.integer(npde), as.integer(kord), as.integer(nint), as.integer(ncc), as.integer(maxder), as.double(X[i]/S[i]), as.double(xbkpt), as.double(work), as.integer(iwork), PACKAGE = "fAsianOptions" ) Price[i] = result[[13]]*S[i] } # Return Value: Price } ################################################################################ # LAPLACE INVERSION: gGemanYor = function(lambda, S = 100, X = 100, Time = 1, r = 0.05, sigma = 0.30, log = FALSE, doplot = FALSE) { # A function written by Diethelm Wuertz # Description: # Calculates function to be Laplace inverted # Arguments: # lambda - complex vector # Notes: # Equation 4.9 with notation as in # Sudler G.F. [1999], "Asian Options: Inverse Laplace # Transform and Martingale Methods Revisited". # FUNCTION: # Settings: x = lambda g = rep(complex(real = 0, imag = 0), length = length(x)) # Calculate for each lambda value from Kummer Function: # Note Kummer function is not vectorized in Indexes ! for ( i in 1:length(x) ) { # Settings: nu = 2*r/(sigma^2) - 1 mu = sqrt(2*lambda[i] + nu^2) q = (sigma^2)*X*Time/(4*S) gamma1 = (mu-nu)/2 gamma2 = (mu+nu)/2 # Convergence Parameters: a = gamma1-2 + 1 b = gamma2+1 + a + 1 z = -1/(2*q) # From Kummer Function [one of ...]: # Use logarithmic Kummer and Gamma Functions to prevent # from numerical overflow! g[i] = kummerM(-z, b-a, b, lnchf = 1) + a*log(-z) + z + cgamma(b-a, log = TRUE) - cgamma(b, log = TRUE) - log (lambda[i]*(lambda[i] - 2 - 2 * nu)) if (!log) g[i] = exp(g[i]) } # Plot function if desired: if (doplot) { if (!is.complex(lambda)) { lam = lambda xlab = "lambda" ylab = "g" } else { lam = Im(lambda) xlab = "Im(lambda)" ylab = "Re(g)" } lambda.min = 4*r/sigma^2 cat("\nmin lambda:", lambda.min, "\n") # Function to be Laplace inverted: print(cbind(lambda, g)) plot(lam, Re(g), type = "l", main = "Laplace Inverse", xlab = xlab, ylab = ylab) lines(lam, 0*Re(g)) # Convergence Indexes: mu = sqrt(2*lambda + nu^2) gamma1 = (mu-nu)/2; a = Re(gamma1-2 + 1) gamma2 = (mu+nu)/2; b = Re(gamma2+1 + a + 1) plot(lam, b, ylim = c(min(c(a,b)), max(c(a,b))), type = "n", xlab = xlab, ylab = "a b", main = "Convergence Indexes") lines(lam, a, col = "red") lines(lam, b, col = "blue") lines(lam, 0*b, type = "l", col = "black") lines(x = lambda.min*c(1,1), y = c( min(c(a,b)), max(c(a,b)) ) ) } # Return Value: g } # ------------------------------------------------------------------------------ GemanYorAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30, doprint = FALSE) { # A function written by Diethelm Wuertz # Description: # Valuate Asian options by Laplace inversion. # FUNCTION: # Parameters: TypeFlag = TypeFlag[1] nu = (2*r)/(sigma^2) - 1 alpha = 1 / sigma^2 * (2*nu + 2) h = sigma^2*Time/4 # Function to be inverted: fx = function(x, S, X, Time, roh, sigma) { nu = (2*roh)/(sigma^2) - 1 alpha = 1 / sigma^2 * (2*nu + 2) h = sigma^2*Time/4 zi = complex(real = 0, imag = x) zc = complex(real = alpha, imag = x) g2 = gGemanYor(lambda = zc, S = S, X = X, Time = Time, r = roh, sigma = sigma, log = TRUE) # Return Value: Re ( exp(zi*h + g2) / (2*pi) ) } # Call: # Integrate stepwise until (hopefully) convergence is reached: q = (sigma^2)*X*Time/(4*S) delta = 10/(2*q) eps = 1.0e-20 if (doprint) { cat("\nDelta:", delta) cat("\nS:", S, "X:", X) cat("\nTime:", Time, "r:", r, "sigma:", sigma, "\n") } i = 1 I = integrate(fx, lower = 0, upper = delta, S = S, X = X, Time=Time, roh = r, sigma = sigma) Price = Increment = exp(-r*Time)*(S/h)*exp(alpha*h)*2*I$value while (abs(Increment)/abs(Price) > eps) { i = i+1 I = integrate(fx, lower = (i-1)*delta, upper = i*delta, S = S, X = X, Time = Time, roh = r, sigma = sigma) Increment = exp(-r*Time)*(S/h)*exp(alpha*h)*2*I$value Price = Price + Increment if (doprint) print(c(i*delta, Price, Increment)) } # Put: # Use Call-Put Parity: if (TypeFlag == "p") { Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X Price = Price - Parity } # Return Value: option = list( price = Price, call = match.call() ) class(option) = "option" option } ################################################################################ # SPECTRAL EXPANSION: gLinetzky = function(x, y, tau, nu, ip = 0) { # A function implemented by Diethelm Wuertz # Description: # Calculates the function to be integrated for the Put Price # in the expression for the spectral representation of # $ P^{(\nu)} (k, \tau) $. Proposition 2 described bu eq. (16) # Note: # Requires Confluent Hypergeometric Functions # Reference: # [L] V. Linetzky, Spectral Expansions for Asian (Average Price) # Options, Preprint, revised Version from October 2002 # Function: result = V = rep(0, length = length(x)) for (i in 1:length(x)) { p = x[i] if (p == 0) { result[i] = 0 } else { z = 1/(2*y) kappa = -(nu+3)/2 mu = complex(real = 0, imag = p/2) # V1 = exp( -(nu^2+p^2)*tau/2 ) * z^kappa * exp(-z/2) logV1 = -(nu^2+p^2)*tau/2 + kappa*log(z) - z/2 # V2 = (abs(cgamma(nu/2+mu)))^2 if (p < 100) { logV2 = log ( (abs(cgamma(nu/2+mu)))^2 ) } else { # Shift by Pi - Take of the proper Phi g = cgamma(nu/2+mu, log = TRUE) r = abs(g) phi = atan(Im(g)/Re(g)) + pi logV2 = 2*r*cos(phi) } # V3 = sinh(pi*p) * p logV3 = pi*p + log(1/2 - exp(-2*pi*p)/2) + log(p) # Whittaker: if (p < 100) { V4 = Re ( whittakerW(z, kappa, mu, ip ) ) } else { # Use: 2 * Re ( (cgamma(-2*mu)/cgamma(1/2-mu-kappa)) * # exp(-z/2) * z^(1/2+mu) * kummerM(z, 1/2+mu-kappa, 1+2*mu) g = log(2) + cgamma(-2*mu, log=TRUE) - cgamma(1/2-mu-kappa, log=TRUE) - z/2 +(1/2+mu)*log(z) + kummerM(z, 1/2+mu-kappa, 1+2*mu, lnchf=1, ip=ip) r = abs(g) # Shift by Pi - Take of the proper Phi phi = atan(Im(g)/Re(g)) + pi logV4 = r*cos(phi) argV4 = cos(r*sin(phi)) } # Collect all terms: if ( p < 100) { result[i] = exp(logV1+logV2+logV3)*V4/(8*pi^2) } else { result[i] = exp(logV1+logV2+logV3+logV4)*argV4/(8*pi^2) } # print(c(p, logV5, logV5a, argV5, argV5a)) } } # Return Value: result } # ------------------------------------------------------------------------------ LinetzkyAsianOption = function(TypeFlag = c("c", "p"), S = 2, X = 2, Time = 1, r = 0.02, sigma = 0.1, table = NA, lower = 0, upper = 100, method = "adaptive", subdivisions = 100, ip = 0, doprint = TRUE, doplot = TRUE,...) { # A function implemented by Diethelm Wuertz # Test for Table: if (is.data.frame(table)) { S = table[,1] X = table[,2] Time = table[,3] r = table[,4] sigma = table[,5] } if (doprint) print(cbind(S, X, Time, r, sigma)) Price = rep(0, length=length(S)) # Settings: tau = sigma^2*Time/4 k = tau*X/S nu = 2*r/sigma^2 - 1 if (doprint) print(cbind(k, tau, nu)) # Parameters: z = 1 / (2*k) kappa = -(nu+3)/2 mu.max = upper/2 if (doprint) print(cbind(z, kappa, mu.max)) # Calculate Spectral Measure: P = P1 = rep(0, times=length(S)) for (i in 1:length(S)) { if (method == "adaptive") { P[i] = integrate(gLinetzky, lower = lower, upper = upper, y = k[i], tau = tau[i], nu = nu[i], ip = ip, subdivisions = subdivisions, rel.tol = .Machine$double.eps^0.25, abs.tol = .Machine$double.eps^0.25, stop.on.error = FALSE)$value } if (method == "trapez") { x = seq(lower, upper, length = subdivisions+1) delta = (upper-lower)/subdivisions F = gLinetzky(x = x, y = k[i], tau[i], nu[i]) # print(c(F[1], F[length(F)], min(F), max(F))) P[i] = ( sum(F)-(F[1]+F[length(F)])/2 ) * delta } if (method == "simpson") { x = seq(lower, upper, length = subdivisions+1) delta = (upper-lower)/subdivisions F = gLinetzky(x = x, y = k[i], tau = tau[i], nu = nu[i], ip = ip) # print(c(F[1], F[length(F)], min(F), max(F))) FF = matrix(F[2:length(F)], byrow = TRUE, ncol = 2) P[i] = (F[1]+4*sum(FF[,1])+2*sum(FF[,2])-F[length(F)]) * delta[i]/3 } # For nu < 0 add: P1[i] = 0 if (nu[i] < 0) { z = 1/(2*k[i]) P1[i] = (2*k[i]*pgamma(abs(nu[i]), z) - pgamma(abs(nu[i])-1, z)) / ( 2 * gamma(abs(nu[i])) ) P[i] = P[i] + P1[i] } # Plot: if (doplot) { x = seq(lower, upper, length = subdivisions) F = gLinetzky(x = x, y = k[i], tau = tau[i], nu = nu[i]) plot(x, F, type = "l") lines(x, 0*x, col = "red") lines(x, F) } } if (doprint) print(cbind(P, P1)) # Derive Call/Put Price: # Put Price: Linetzky = exp(-r*Time) * (S/tau) * P # Call Price: if (TypeFlag == "c") { Linetzky = Linetzky + S*(1-exp(-r*Time))/(r*Time) - X*exp(-r*Time) } # Return Value: option = list( price = Linetzky, call = match.call() ) class(option) = "option" option } ################################################################################ # ASIAN BOUNDS: # Note: # We have not implemented the formula for the upper bound derived # by Rogers and Shi. Thompson's upper bound formula is much more # precise and therefore we have concentrated ourself on their # approach. BoundsOnAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30, table = NA, method = c("CT", "RST", "T")) { # A function implemented by Diethelm Wuertz # Description: # Calculates Bounds on Asian Option Prices # CT - Curran-Thompson Lower Bound # RST - Roger-Shi-Thompson Lower Bound # T - Thompson Upper Bound # FUNCTION: # Set Default Method, if no other is selected: TypeFlag = TypeFlag[1] if (length(method) == 3) method = "T" # Test for Table: if (is.data.frame(table)) { S = table[,1] X = table[,2] Time = table[,3] r = table[,4] sigma = table[,5] } Price = rep(NA, length = length(S)) # Curran-Thompson Lower Bound: if (method == "CT") { for ( i in 1:length(S) ) { Price[i] = CurranThompsonAsianOption(TypeFlag=TypeFlag, S=S[i], X=X[i], Time[i], r=r[i], sigma[i])$price } } # Roger-Shi-Thompson Lower Bound: if (method == "RST") { for ( i in 1:length(S) ) { Price[i] = RogerShiThompsonAsianOption(TypeFlag=TypeFlag, S = S[i], X = X[i], Time[i], r = r[i], sigma[i])$price } } # Thompson Upper Bound: if (method == "T") { for ( i in 1:length(S) ) { Price[i] = ThompsonAsianOption(TypeFlag = TypeFlag, S = S[i], X = X[i], Time[i], r = r[i], sigma[i])$price } } # Return Value: option = list( price = Price, call = match.call() ) class(option) = "option" option } # ------------------------------------------------------------------------------ CurranThompsonAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Description: # Calculates "lower bound" for Asian Call Option from # Thompson's formula describing the continuous limit # of Curran's approximation. # Note: # Rescale sigma: # Note the formula of Thompson work for Time=1 only! # Thus the easiest way to cover times to maturity different # from unity can be achieved by scale the volatility and # interest rate! - Just do it # FUNCTION: # Settings: TypeFlag = TypeFlag[1] sigma = sigma*sqrt(Time) r = r*Time Time = 1 # Settings: alpha = r - 0.5*sigma^2 # Solve for gamma star: f1 = function(x, S, X, alpha, sigma) { exp( 3 * ( log(X/S)-alpha/2 ) * x *(1-x/2) + alpha * x + 0.5 * sigma^2 * (x-3*x^2*(1-x/2)^2) ) } # Integrate: gs = integrate(f1, lower = 0, upper = 1, S = S, X = X, alpha = alpha, sigma = sigma, subdivisions = 1000, rel.tol = .Machine$double.eps^0.5, abs.tol = .Machine$double.eps^0.5)$value # Final Calculate: g.star = ( log(2*X/S-gs) - alpha/2 ) / sigma # Solve for lower bound: f = function(x, g.star, alpha, sigma) { time = x arg = (-g.star + sigma*time*(1-time/2))*sqrt(3) f = exp( (alpha+sigma^2/2)*time ) * pnorm(arg) f } # Integrate: value = integrate(f, lower=0, upper=1, g.star=g.star, alpha=alpha, sigma=sigma, subdivisions=1000, rel.tol=.Machine$double.eps^0.5, abs.tol=.Machine$double.eps^0.5)$value # Call Price: CurranThompson = exp(-r*Time) * ( S*value - X*pnorm(-g.star*sqrt(3)) ) # Put Price: if (TypeFlag == "p") { Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X CurranThompson = CurranThompson - Parity } # Return Value: option = list( price = CurranThompson, call = match.call() ) class(option) = "option" option } # ------------------------------------------------------------------------------ RogerShiThompsonAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Description: # Calculates "lower bound" for Asian Call Option from # Thompson's formula. Thompson's result is the same # as can be obtained from Roger and Shi's formula. # However, Thompson's formula is numerically more # efficient since it requires a single integration # only, whereas Roger and Shi's formula requires # double integration. # Note: # Rescale sigma: # Note the formula of Thompson work for Time=1 only! # Thus the easiest way to cover times to maturity different # from unity can be achieved by scale the volatility and # interest rate! - Just do it # FUNCTION: # Settings: TypeFlag = TypeFlag[1] sigma = sigma*sqrt(Time) r = r*Time Time = 1 # Function from which to calculate gamma star: gamma.star = function(S, X, r, sigma, lower = -99, upper = 99) { func = function(gamma, S, X, r, sigma) { f = function(x, gamma, roh, sigma) { time = x alpha = roh - sigma*sigma/2 f = exp( 3 * gamma * sigma * time * (1-time/2) + alpha * time + sigma * sigma * (time-3*time^2*(1-time/2)^2)/2 ) f } # Integrate: integrate(f, lower = 0, upper = 1, gamma = gamma, roh = r, sigma = sigma, subdivisions = 1000, rel.tol = .Machine$double.eps^0.5, abs.tol = .Machine$double.eps^0.5)$value - X/S } # Find Root Value: uniroot(func, lower = lower, upper = upper, S = S, X = X, r = r, sigma = sigma)$root } g.star = gamma.star(S, X, r, sigma) # Function to be integrated: f = function(x, g.star, roh, sigma) { time = x alpha = roh - sigma*sigma/2 arg = (-g.star + sigma*time*(1-time/2))*sqrt(3) f = exp( (alpha+sigma^2/2)*time ) * pnorm(arg) f } # Integrate: value = integrate(f, lower=0, upper=1, g.star=g.star, roh=r, sigma=sigma, subdivisions=1000, rel.tol=.Machine$double.eps^0.5, abs.tol=.Machine$double.eps^0.5)$value # Call Price: RogerShiThompson = exp(-r*Time) * ( S*value - X*pnorm(-g.star*sqrt(3)) ) # Put Price: if (TypeFlag == "p") { Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X RogerShiThompson = RogerShiThompson - Parity } # Return Value: option = list( price = RogerShiThompson, call = match.call() ) class(option) = "option" option } # ------------------------------------------------------------------------------ ThompsonAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Description: # Calculates "upper bound" for Asian Call Option from # Thompson's formula. # Note: # Rescale sigma: # Note the formula of Thompson work for Time=1 only! # Thus the easiest way to cover times to maturity different # from unity can be achieved by scale the volatility and # interest rate! - Just do it # FUNCTION: # Settings: TypeFlag = TypeFlag[1] sigma = sigma*sqrt(Time) r = r*Time Time = 1 # Internal Functions: sqrtvt = function(x, S, X, alpha, sigma) { t = x ct = S*exp(alpha*t)*sigma - X*sigma vt = ct^2*t + 2*(X*sigma)*ct*t*(1-t/2) + (X*sigma)^2/3 sqrt(vt) } fmu = function(t, S, X, r, sigma) { alpha = r -sigma^2/2 gint = integrate(sqrtvt, lower=0, upper=1, S=S, X=X, alpha=alpha, sigma=sigma, subdivisions=1000, rel.tol=.Machine$double.eps^0.5, abs.tol=.Machine$double.eps^0.5)$value gamma = (X - S * (exp(alpha)-1) / alpha ) / gint mu = (S*exp(alpha*t) + gamma*sqrtvt(x=t, S=S, X=X, alpha=alpha, sigma=sigma) ) / X mu } fatx = function(t, x, S, X, r, sigma) { alpha = r - sigma^2/2 mu = fmu(t=t, S=S, X=X, r=r, sigma=sigma) atx = S*exp(sigma*x+alpha*t) - X*(mu + sigma*x) + X*sigma*(1-t/2)*x atx } fbtx = function(t, x, S, X, r, sigma) { alpha = r - sigma^2/2 btx = X*sigma*sqrt(1/3 -t*(1-t/2)^2) btx } fw = function(x, v, S2, X2, r2, sigma2) { w = x atx = fatx(t=v^2, x=w*v, S=S2, X=X2, r=r2, sigma=sigma2) btx = fbtx(t=v^2, x=w*v, S=S2, X=X2, r=r2, sigma=sigma2) 2 * v * dnorm(w) * (atx*pnorm(atx/btx) + btx*dnorm(atx/btx)) } fv = function(x, S1, X1, r1, sigma1) { fv = rep(0, length=length(x)) for (i in 1:length(x)) fv[i] = integrate(fw, lower=-20, upper=20, v=x[i], S2=S1, X2=X1, r2=r1, sigma2=sigma1, subdivisions=1000, rel.tol=.Machine$double.eps^0.5, abs.tol=.Machine$double.eps^0.5)$value exp(-r)*fv } # Integrate - Call Price: Thompson = integrate(fv, lower=0, upper=1, S1=S, X1=X, r1=r, sigma1=sigma, subdivisions=1000, rel.tol=.Machine$double.eps^0.5, abs.tol=.Machine$double.eps^0.5)$value # Put Price: if (TypeFlag == "p") { Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X Thompson = Thompson - Parity } # Return Value: option = list( price = Thompson, call = match.call() ) class(option) = "option" option } # ------------------------------------------------------------------------------ TolmatzAsianOption = function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz # Description: # Calculates "lower bound" for Asian Call Option from # the asymptotic behavior derived by Tolmatz # FUNCTION: TypeFlag = TypeFlag[1] Tolmatz = NA # Put Price: if (TypeFlag == "p") { Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X Tolmatz = Tolmatz - Parity } # Return Value: option = list( price = Tolmatz, call = match.call() ) class(option) = "option" option } ################################################################################ # SYMMETRY AND EQUIVALENCE RELATIONS: CallPutParityAsianOption = function(TypeFlag = "p", Price = 8.828759, S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.3, table = NA) { # A function implemented by Diethelm Wuertz # Test for Table: if (is.data.frame(table)) { S = table[,1] X = table[,2] Time = table[,3] r = table[,4] sigma = table[,5] } # Call from Put: if (TypeFlag == "c") { # print(Price) Parity = S*(1-exp(-r*Time))/(r*Time) - X*exp(-r*Time) # print(Parity) result = Price + Parity } # Put from Call: if (TypeFlag == "p") { Parity = S*(1-exp(-r*Time))/(r*Time) - X*exp(-r*Time) result = Price - Parity } # Return Value: result } # ------------------------------------------------------------------------------ WithDividendsAsianOption = function(TypeFlag = "c", Dividends = 0.45, S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.3, calculator = MomentMatchedAsianOption, method = "LN") { # A function implemented by Diethelm Wuertz # Add Dividends: q = Dividends = 0.05 r.q = r - q X.q = X * exp(-q*Time) S.q<- S * exp(-q*Time) # Call Price: if (TypeFlag == "c") Price = calculator(TypeFlag = TypeFlag, S = S.q, X = X.q, Time = Time, r = r.q, sigma = sigma, method = method)$price # Put Price: if (TypeFlag == "p" ) Price = NA # Return Value: option = list( price = Price, call = match.call() ) class(option) = "option" option } ################################################################################ # TABULATED RESULTS: FuMadanWangTable = function() { # A function implemented by Diethelm Wuertz # Description: # Display Fu-Madan-Wang's results from Table 1 # Source: # Settings: X = rep(c(90,95,100,105,110), times = 6) S = 100*rep(1, times = length(X)) sigma = rep(0.20, length = length(X)) r = rep(0.09, length = length(X)) # Call Prices: CallEulerFMW = c( 11.5293, 7.2131, 3.8087, 1.6465, 0.5761, 11.9247, 7.7249, 4.3696, 2.1175, 0.8734, 13.8372, 9.9998, 6.7801, 4.2982, 2.5473, 17.1212, 13.6763, 10.6319, 8.0436, 5.9267, 19.8398, 16.6740, 13.7974, 11.2447, 9.0316, 24.0861, 21.3774, 18.8399, 16.4917, 14.3442) CPUEulerFMW = c( 53,44,42,45,38, 34,33,33,32,32, 26,26,26,25,25, 24,23,23,23,22, 21,21,21,20,20, 19,22,19,21,21) CallPostWidderFMW = c( 11.5176, 7.1981, 3.8196, 1.6623, 0.5728, 11.9241, 7.7185, 4.3759, 2.1290, 0.8753, 13.8439, 10.0029, 6.7823, 4.3010, 2.5450, 17.1297, 13.6830, 10.6370, 8.0474, 5.9295, 19.8495, 16.6822, 13.8042, 11.2502, 9.0360, 24.0861, 21.3774, 18.8399, 16.4917, 14.3442) CPUPostWidderFMW = c( 640,631,628,623,625, 613,591,601,585,595, 518,513,514,503,503, 475,473,473,469,474, 468,467,467,462,465, 453,442,449,441,453) CallZhang = c( NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 13.8314996, 9.99566567, 6.7773481, 4.2964626, 2.5462209, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA) # Return Value: data.frame(cbind( S, X, r, sigma, CallEulerFMW, CPUEulerFMW, CallPostWidderFMW, CPUPostWidderFMW, CallZhang)) } # ------------------------------------------------------------------------------ FusaiTaglianiTable = function() { # A function implemented by Diethelm Wuertz # Description: # Display Fusai and Tagliani's results from Table 6a - 6c # Source: # G. Fusai and A. Tagliani [2002] # An Accurate Valuation of Asian Options Using Moments # Settings: S = 100*rep(1, times = 45) X = rep(rep(c(90, 95, 100, 105, 110), times = 3), times = 3) Time = rep(1, times = 45) r = rep(sort(rep(c(0.05, 0.09, 0.15), times = 5)), times = 3) sigma = rep(0.10, times = 15); sigma = c(sigma, 3*sigma, 5*sigma) # Moment Matched Call Prices: CallLNa = c( 11.95333, 7.41517, 3.64748, 1.30684, 0.32399, 13.38629, 8.91721, 4.92310, 2.07045, 0.62338, 15.39906, 11.12196, 7.03485, 3.61869, 1.41080) CallLNb = c( 14.03812, 10.74879, 7.99251, 5.77417, 4.05724, 15.06704, 11.73287, 8.88576, 6.54628, 4.69511, 16.59082, 13.22661, 10.27814, 7.78408, 5.74800) CallLNc = c( 17.67918, 14.91574, 12.48954, 10.38535, 8.58075, 18.43698, 15.66486, 13.21198, 11.06751, 9.21323, 19.55391, 16.78229, 14.30234, 12.10904, 10.18997) CallFTLN = c(CallLNa, CallLNb, CallLNc) # Gram-Charlier Call Prices: CallFTa = c( 11.95127, 7.40754, 3.64091, 1.31097, 0.33156, 13.38535, 8.91185, 4.91459, 2.07002, 0.63072, 15.39885, 11.11966, 7.02746, 3.61216, 1.41384) CallFTb = c( 13.89562, 10.63393, 7.92948, 5.76852, 4.09909, 14.92543, 11.60605, 8.80190, 6.51749, 4.71737, 16.45856, 13.09056, 10.16897, 7.72184, 5.73774) CallFTc = c( 17.00382, 14.46257, 12.25646, 10.34604, 8.69620, 17.69986, 15.13557, 12.89937, 10.95396, 9.26553, 18.73925, 16.14761, 13.87252, 11.88030, 10.13926) CallFTLNGC = c(CallFTa, CallFTb, CallFTc) # Return Value: data.frame(S, X, Time, r, sigma, CallFTLN, CallFTLNGC) } # ------------------------------------------------------------------------------ GemanTable = function() { # A function implemented by Diethelm Wuertz # Description: # Display Geman's Table # Source: # H. Geman [], # Functionals of Brownian Motion in Path Dependent # Option Valuation # Settings: S = c(1.9, 2.0, 2.1, 2.0, 2.0, 2.0) X = rep(2, times = 6) Time = c(1, 1, 1, 1, 2, 2) r = c(5, 5, 5, 2, 1.25, 5)/100 sigma = c(5, 5, 5, 1, 2.5, 5)/10 # Call Prices: CallGY = c( 0.195, 0.248, 0.308, 0.058, 0.1772, 0.352) CallMC = c( 0.191, 0.248, 0.306, 0.056, 0.1771, 0.347) # Return Value: data.frame(cbind(S, X, Time, r, sigma, CallGY, CallMC)) } # ------------------------------------------------------------------------------ LinetzkyTable = function() { # A function implemented by Diethelm Wuertz # Description: # Display Linetzky's Table 3 # Source: # V. Linetzky [2002] # Spectral Expansions for Asian (Average Price) Options # Settings: S = c(2.00, 2.00, 2.00, 1.90, 2.00, 2.10, 2.00) X = c(2.00, 2.00, 2.00, 2.00, 2.00, 2.00, 2.00) Time = c(1.00, 1.00, 2.00, 1.00, 1.00, 1.00, 2.00) r = c(0.02, 0.18, 0.0125, 0.05, 0.05, 0.05, 0.05) sigma = c(10.0, 30.0, 25.0, 50.0, 50.0, 50.0, 50.0) / 100 # Call Prices: CallEE = c( 0.0559860415, 0.2183875466, 0.1722687410, 0.1931737903, 0.2464156905, 0.3062203648, 0.3500952190) CallSLT = c( 0.055986, 0.218388, 0.172269, 0.193174, 0.246416, 0.306220, 0.350095) CallMC = c( 0.05602, 0.2185, 0.1725, 0.1933, 0.2465, 0.3064, 0.3503) CallTLB = c( 0.055985, 0.218366, 0.172226, 0.193060, 0.246298, 0.306094, 0.349779) CallTUB = c( 0.055989, 0.218473, 0.172451, 0.193799, 0.247054, 0.306904, 0.352556) # Return Value: data.frame(S, X, Time, r, sigma, CallEE, CallSLT, CallMC, CallTLB, CallTUB) } # ------------------------------------------------------------------------------ ZhangTable = function() { # A function implemented by Diethelm Wuertz # Description: # Display Zhangs's results from Table 1 # Source: # J.E. Zhang [2002], # A semi-analytical method for pricing and hedging # continuously sampled arithmetic average rate options # Settings: S = 100*rep(1, times = 36) X = c(rep(c(95, 100, 105), times = 3), rep(c(90, 100, 110), times = 9)) Time = rep(1, times = 36) r = rep(c(5, 5, 5, 9, 9, 9, 15, 15, 15)/100, times = 4) sigma = sort(rep(c(0.05, 0.10, 0.20, 0.30), times = 9)) # Call Prices: CallZ = c( 7.1777275, 2.7161745, 0.3372614, 8.8088302, 4.3082350, 0.9583841, 11.0940944, 6.7943550, 2.7444531, 11.9510927, 3.6413864, 0.3312030, 13.3851974, 4.9151167, 0.6302713, 15.3987687, 7.0277081, 1.4136149, 12.5959916, 5.7630881, 1.9898945, 13.8314996, 6.7773481, 2.5462209, 15.6417575, 8.4088330, 3.5556100, 13.9538233, 7.9456288, 4.0717942, 14.9839595, 8.8287588, 4.6967089, 16.5129113, 10.2098305, 5.7301225) # Return Value: data.frame(cbind(S, X, Time, r, sigma, CallZ)) } # ------------------------------------------------------------------------------ ZhangShortTable = function() { # A function implemented by Diethelm Wuertz # Description: # Display Zhangs's short-tenor (ST) results from Table 7 # Source: # J.E. Zhang [2002], # A semi-analytical method for pricing and hedging # continuously sampled arithmetic average rate options # Settings: S = 100*rep(1, times = 18) X = c(rep(c(95,100,105), times = 6)) Time = rep(0.08, times = 18) r = rep(0.09, times = 18) sigma = sort(rep(c(0.05, 0.10, 0.20, 0.30, 0.40, 0.50), times = 3)) # Call Prices: CallPM = c( 5.3224059, 0.5343219, 0.0000000, 5.3225549, 0.8431357, 0.0014921, 5.3816262, 1.4835707, 0.1292560, 5.6304554, 2.1291126, 0.4880727, 6.0222502, 2.7758194, 0.9753549, 6.4947818, 3.4228569, 1.5274729) CallCT = c( 5.3224059, 0.5343040, 0.0000000, 5.3225550, 0.8431302, 0.0014924, 5.3816340, 1.4835272, 0.1292529, 5.6304480, 2.1289662, 0.4879999, 6.0221516, 2.7754740, 0.9750983, 6.4944850, 3.4221863, 1.5268836) # Return Value: data.frame(cbind(S, X, Time, r, sigma, CallPM, CallCT)) } # ------------------------------------------------------------------------------ ZhangLongTable = function() { # A function implemented by Diethelm Wuertz # Description: # Display Zhangs's long-tenor (LT) results from Table 7 # Source: # J.E. Zhang [2002], # A semi-analytical method for pricing and hedging # continuously sampled arithmetic average rate options # Settings: S = 100*rep(1, times = 18) X = c(rep(c(95,100,105), times = 6)) Time = rep(3, times = 18) r = rep(0.09, times = 18) sigma = sort(rep(c(0.05, 0.10, 0.20, 0.30, 0.40, 0.50), times = 3)) # Call Prices: CallPM = c( 15.11626, 11.30359, 7.55327, 15.21331, 11.63716, 8.39113, 16.63363, 13.76592, 11.22170, 19.01304, 16.58331, 14.39723, 21.70848, 19.57038, 17.62156, 24.47434, 22.55837, 20.79507) CallCT = c( 15.11626, 11.30361, 7.55328, 15.21376, 11.63752, 8.39084, 16.63611, 13.76547, 11.21783, 19.01856, 16.58101, 14.38708, 21.72991, 19.57593, 17.61228, 24.54788, 22.60645, 20.81811) # Return Value: data.frame(cbind(S, X, Time, r, sigma, CallPM, CallCT)) } ################################################################################