diff -Nru fasianoptions-2100.76/ChangeLog fasianoptions-3000.78/ChangeLog --- fasianoptions-2100.76/ChangeLog 2009-09-28 15:21:14.000000000 +0000 +++ fasianoptions-3000.78/ChangeLog 2013-02-07 21:34:56.000000000 +0000 @@ -0,0 +1,52 @@ +2013-02-07 chalabi + + * R/zzz.R: Removed deprecated .First.lib() + * R/GammaFunctions.R: Changed deprecated is.real() by is.double() + * DESCRIPTION: Updated maintainer field + +2012-03-21 chalabi + + * ChangeLog, DESCRIPTION: updated DESCRIPTION and ChangeLog + * src/EBMAsianOptions.f: replaced 'write' calls by R Fortran print + API + +2012-03-20 chalabi + + * src/EBMAsianOptions.f: fixed gfortran warnings in Fortran + routines + * DESCRIPTION: updated DESC file + +2012-03-19 chalabi + + * R/EBMAsianOptions.R, R/GammaFunctions.R, + R/HypergeometricFunctions.R: removed partial argument names + * src/EBMAsianOptions.f: removed part of the WRITE() statements in + Fortran code + * NAMESPACE: added NAMESPACE + +2011-09-23 mmaechler + + * DESCRIPTION: remove deprecated "LazyLoad" entry + +2010-07-23 chalabi + + * inst/DocCopying.pdf: removed DocCopying.pdf license is already + specified in DESCRIPTION file + +2009-10-01 chalabi + + * DESCRIPTION: updated version number + +2009-09-29 chalabi + + * ChangeLog, DESCRIPTION: updated DESC and ChangeLog + +2009-04-02 chalabi + + * DESCRIPTION: more explicit depends and suggests field in DESC + file. + +2009-04-01 chalabi + + * DESCRIPTION: updated DESC file + diff -Nru fasianoptions-2100.76/DESCRIPTION fasianoptions-3000.78/DESCRIPTION --- fasianoptions-2100.76/DESCRIPTION 2009-09-30 19:40:52.000000000 +0000 +++ fasianoptions-3000.78/DESCRIPTION 2013-02-08 16:03:21.000000000 +0000 @@ -1,21 +1,20 @@ Package: fAsianOptions -Version: 2100.76 -Revision: 4057 -Date: 2009-09-28 +Version: 3000.78 +Revision: 5443 +Date: 2013-02-07 Title: EBM and Asian Option Valuation Author: Diethelm Wuertz and many others, see the SOURCE file Depends: R (>= 2.4.0), timeDate, timeSeries, fBasics, fOptions Suggests: RUnit -Maintainer: Rmetrics Core Team +Maintainer: Yohan Chalabi Description: Environment for teaching "Financial Engineering and Computational Finance" -NOTE: SEVERAL PARTS ARE STILL PRELIMINARY AND MAY BE CHANGED IN THE - FUTURE. THIS TYPICALLY INCLUDES FUNCTION AND ARGUMENT NAMES, AS - WELL AS DEFAULTS FOR ARGUMENTS AND RETURN VALUES. -LazyLoad: yes +Note: Several parts are still preliminary and may be changed in the + future. this typically includes function and argument names, as + well as defaults for arguments and return values. LazyData: yes License: GPL (>= 2) URL: http://www.rmetrics.org -Packaged: 2009-09-28 15:21:34 UTC; yankee +Packaged: 2013-02-08 15:12:54 UTC; yankee Repository: CRAN -Date/Publication: 2009-09-30 19:40:52 +Date/Publication: 2013-02-08 17:03:21 diff -Nru fasianoptions-2100.76/MD5 fasianoptions-3000.78/MD5 --- fasianoptions-2100.76/MD5 1970-01-01 00:00:00.000000000 +0000 +++ fasianoptions-3000.78/MD5 2013-02-08 16:03:21.000000000 +0000 @@ -0,0 +1,27 @@ +32bf4eb915a4a8d904030771c17c0dc6 *ChangeLog +c0daa94321e843e245dd3263d17b2387 *DESCRIPTION +5d0d38953901cdabfa1c56977da4fdf1 *NAMESPACE +be7df6130ace1097ba6e8ac0971a7089 *R/BesselFunctions.R +b24a5a867d0ecc11c159b1b58cc40393 *R/EBMAsianOptions.R +9882507e5c37d8a890edcaf01bf6f696 *R/EBMDistribution.R +e7efd48788e0248853efb338c3fa578e *R/GammaFunctions.R +79480e99627d3c023f20c88e367debbf *R/HypergeometricFunctions.R +3c77bef693ed7c293f2a3319d1edbdba *R/zzz.R +6042b9c5e5bec3ecc1b6959cd2858b64 *inst/COPYRIGHT.html +b208277c86926092266ea1ba1d2d9e0b *inst/unitTests/Makefile +a516e208a0f1ef3c9f8e2f5e9b36290a *inst/unitTests/runTests.R +797b2a14d90a119ec5660bc95ef33610 *inst/unitTests/runit.BesselFunctions.R +8f4de115f88501c9f00d45ca58fa9fe6 *inst/unitTests/runit.EBMAsianOptions.R +fa7baf4aa65c517a5c53548453af02a8 *inst/unitTests/runit.EBMDistribution.R +96610bc14fb29e4680144e0efe31ef33 *inst/unitTests/runit.GammaFunctions.R +197f80aaca6c6221faaabb88d40b083b *inst/unitTests/runit.HypergeometricFunctions.R +c0d59ee25dc707e515d2c4953bf1e630 *man/BesselFunctions.Rd +bb017f4ecab6bab8d23b1fec1d06c159 *man/EBMAsianOptions.Rd +b898124e76be416234c45328f391d2ef *man/EBMDistribution.Rd +b206863027bdd6b6f48e264f63352785 *man/GammaFunctions.Rd +b0397e3bbfb223909451b1ed199b6894 *man/HypergeometricFunctions.Rd +a741df3291ffa772d712e7660ac0c66f *src/EBMAsianOptions.f +00b5c600bd679e0127d91a08cddcbe72 *src/GammaFunctions.f +02b7119a53ace717d00d3538d8ffdd41 *src/HypergeometricFunctions.f +1d848c9e918f0eafbb64d4aa0e593e78 *src/Makevars +ca566e590ec30abd0718c5375e1a446f *tests/doRUnit.R diff -Nru fasianoptions-2100.76/NAMESPACE fasianoptions-3000.78/NAMESPACE --- fasianoptions-2100.76/NAMESPACE 1970-01-01 00:00:00.000000000 +0000 +++ fasianoptions-3000.78/NAMESPACE 2012-03-19 10:06:08.000000000 +0000 @@ -0,0 +1,110 @@ + +################################################ +## import name space +################################################ + +import("timeDate") +import("timeSeries") +import("fBasics") +import("fOptions") + +################################################ +## useDynLib +################################################ + +useDynLib("fAsianOptions") + +################################################ +## S4 classes +################################################ + + +################################################ +## S3 classes +################################################ + + +################################################ +## functions +################################################ + +export( + ".AbrahamsonAsianOptionMoments", + ".Bessel.ENVJ", + ".Bessel.MSTA1", + ".Bessel.MSTA2", + ".Bessel01", + ".BesselN", + ".DufresneAsianOptionMoments", + ".DufresneMoments", + ".GramCharlierAsianDensity", + ".LevyTurnbullWakemanAsianDensity", + ".LevyTurnbullWakemanAsianOption", + ".MilevskyPosnerAsianDensity", + ".MilevskyPosnerAsianOption", + ".PosnerMilevskyAsianDensity", + ".PosnerMilevskyAsianOption", + ".TolmatzAsianOptionMoments", + ".TurnbullWakemanAsianOptionMoments", + ".gxtEBM", + ".gxtuEBM", + ".gxuEBM", + ".psiEBM", + ".thetaEBM", + "AsianOptionMoments", + "BesselDI", + "BesselDK", + "BesselI", + "BesselK", + "BoundsOnAsianOption", + "CallPutParityAsianOption", + "CurranThompsonAsianOption", + "FuMadanWangTable", + "FusaiTaglianiTable", + "GemanTable", + "GemanYorAsianOption", + "GramCharlierAsianOption", + "LinetzkyAsianOption", + "LinetzkyTable", + "MomentMatchedAsianDensity", + "MomentMatchedAsianOption", + "Pochhammer", + "Psi", + "RogerShiThompsonAsianOption", + "ThompsonAsianOption", + "TolmatzAsianOption", + "VecerAsianOption", + "WithDividendsAsianOption", + "ZhangApproximateAsianOption", + "ZhangAsianOption", + "ZhangLongTable", + "ZhangShortTable", + "ZhangTable", + "cgamma", + "d2EBM", + "dEBM", + "dasymEBM", + "derivative", + "dgam", + "djohnson", + "dlognorm", + "drgam", + "erf", + "gGemanYor", + "gLinetzky", + "hermiteH", + "igamma", + "kummerM", + "kummerU", + "masian", + "mjohnson", + "mlognorm", + "mnorm", + "mrgam", + "pEBM", + "pgam", + "pjohnson", + "plognorm", + "prgam", + "whittakerM", + "whittakerW" ) diff -Nru fasianoptions-2100.76/R/EBMAsianOptions.R fasianoptions-3000.78/R/EBMAsianOptions.R --- fasianoptions-2100.76/R/EBMAsianOptions.R 2009-09-02 20:27:40.000000000 +0000 +++ fasianoptions-3000.78/R/EBMAsianOptions.R 2012-03-19 10:08:22.000000000 +0000 @@ -15,7 +15,7 @@ # MA 02111-1307 USA # Copyrights (C) -# for this R-port: +# for this R-port: # 1999 - 2004, Diethelm Wuertz, GPL # Diethelm Wuertz # info@rmetrics.org @@ -36,7 +36,7 @@ # MomentMatchedAsianDensity Valuate moment matched option densities # .LevyTurnbullWakemanAsianDensity Log-Normal Approximation # .MilevskyPosnerAsianDensity Reciprocal-Gamma Approximation -# .PosnerMilevskyAsianDensity Johnson Type I Approximation +# .PosnerMilevskyAsianDensity Johnson Type I Approximation # GRAM CHARLIER SERIES EXPANSION: DESCRIPTION: # GramCharlierAsianOption Calculate Gram-Charlier option prices # .GramCharlierAsianDensity NA @@ -58,7 +58,7 @@ # PDEAsianOption PDE Asian Option Pricing # .ZhangAsianOption Asian option price by Zhang's 1D PDE # ZhangApproximateAsianOption -# .VecerAsianOption Asian option price by Vecer's 1D PDE +# .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 @@ -68,8 +68,8 @@ # 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 +# 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 @@ -88,8 +88,8 @@ # MOMENT MATCHING: -MomentMatchedAsianOption = -function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, +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 @@ -98,132 +98,132 @@ # 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] + 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]) ) ) + 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], + moments = masian(Time = Time[i], r = r[i], sigma = sigma[i])$rawMoments - moments = moments / Time[i]^(1:4) + 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) - + 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]) + 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) + 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 ) - + 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 + Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X + Price = Price - Parity } - + # Return Value: option = list( - price = Price, + price = Price, call = match.call() ) class(option) = "option" option -} +} # ------------------------------------------------------------------------------ -.LevyTurnbullWakemanAsianOption = -function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, +.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") + 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, +.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") + 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, +.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") + MomentMatchedAsianOption(TypeFlag[1], S, X, Time, r, sigma, method = "JI") } # ------------------------------------------------------------------------------ -MomentMatchedAsianDensity = +MomentMatchedAsianDensity = function(x, Time = 1, r = 0.09, sigma = 0.30, method = c("LN", "RG", "JI")) { # A function implemented by Diethelm Wuertz @@ -232,46 +232,46 @@ # 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 = 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) + 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) + 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 = masian(Time = Time, r = r, sigma = sigma)$centralMoments cmoments = cmoments / Time^(1:4) - mu = cmoments[1] - varsigma = sqrt(cmoments[2]) + 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) + 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 } @@ -280,7 +280,7 @@ # ------------------------------------------------------------------------------ -.LevyTurnbullWakemanAsianDensity = +.LevyTurnbullWakemanAsianDensity = function(x, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz @@ -292,7 +292,7 @@ # ------------------------------------------------------------------------------ -.MilevskyPosnerAsianDensity = +.MilevskyPosnerAsianDensity = function(x, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz @@ -304,7 +304,7 @@ # ------------------------------------------------------------------------------ -.PosnerMilevskyAsianDensity = +.PosnerMilevskyAsianDensity = function(x, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz @@ -317,20 +317,20 @@ # GRAM CHARLIER: -GramCharlierAsianOption = -function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, +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 - +{ # A function implemented by Diethelm Wuertz + # Description: - # Calculate arithmetic Asian options price using + # 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] @@ -341,58 +341,58 @@ X = table[,2] Time = table[,3] r = table[,4] - sigma = table[,5] + sigma = table[,5] } - + # Calculate Price: - Price = MomentMatchedAsianOption("c", S = S, X = X, Time = Time, r = r, + 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], + 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 + 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: + + # 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], + 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 + 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]) + 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) + d = sign(eta) c = d*mu - exp( (0.5/b-a)/b ) asian.cm = cmoments johnson.cm = cmoments @@ -400,33 +400,33 @@ 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) + 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 + 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) - + 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 + Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X + Price = Price - Parity } - + # Return Value: option = list( - price = Price, + price = Price, call = match.call() ) class(option) = "option" option -} +} -.GramCharlierAsianDensity = +.GramCharlierAsianDensity = function(Time = 1, r = 0.09, sigma = 0.30, method = c("LN", "RG", "JI")) -{ # A function ported by Diethelm Wuertz +{ # A function ported by Diethelm Wuertz # Return Value: NA @@ -437,8 +437,8 @@ # STATE SPACE MOMENTS: -AsianOptionMoments = -function(M = 4, Time = 1, r = 0.045, sigma = 0.30, log = FALSE, +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 @@ -449,66 +449,66 @@ # TW - First 2 Moments from Turnbull-Wakeman # T - Asymptotic Behavior after Tolmatz - # FUNCTION: - + # FUNCTION: + # Settings: method = method[1] result = NA - + # Abrahamson Formula: - if (method == "A") result = - .AbrahamsonAsianOptionMoments(M = M, Time = Time, r = r, + 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, + 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, + 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, + if (method == "TW") result = + .TurnbullWakemanAsianOptionMoments(M = M, Time = Time, r = r, sigma = sigma) - + # Return Value: - result + result } # ------------------------------------------------------------------------------ -.DufresneAsianOptionMoments = +.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) { + 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 - } + d + } moments = rep(0, length=M) - for (n in 1:M) { + for (n in 1:M) { moments[n] = 0 - for (j in 0:n) moments[n] = moments[n] + + 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 + moments } - + # Calculate for: tau = sigma^2*Time/4 nu = 2*r/sigma^2-1 @@ -517,28 +517,28 @@ (4/sigma^2)^(1:M) * moments(M, tau, nu) } - + # ------------------------------------------------------------------------------ -.AbrahamsonAsianOptionMoments = +.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 = 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]) + if (j!= m) d[m] = d[m]*(c[m]-c[j]) } d[m] = exp(c[m]*tau) / d[m] } @@ -546,42 +546,42 @@ } 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 = +.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 + # by Turnbull and Wakeman. It can serve as a test for # other implementations. - + # Note: # Maximum M is 2! - - # FUNCTION: - + + # FUNCTION: + # Moments: - moments = rep(NA, times = M) + 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) + + 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 } @@ -590,7 +590,7 @@ # ------------------------------------------------------------------------------ -.TolmatzAsianOptionMoments = +.TolmatzAsianOptionMoments = function (M = 100, Time = 1, r = 0.045, sigma = 0.30, log = FALSE) { # A function implemented by Diethelm Wuertz @@ -598,25 +598,25 @@ # 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 = 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 + 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 @@ -643,48 +643,48 @@ # 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) +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 + # 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] + 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 + + # 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)) * @@ -692,7 +692,7 @@ 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", @@ -721,100 +721,100 @@ as.double(xbkpt), as.double(work), as.integer(iwork), - PACKAGE = "fAsianOptions" + PACKAGE = "fAsianOptions" ) Price[i] = result[[13]]*S[i] } - + # ? - Price = Price + - ZhangApproximateAsianOption(TypeFlag, S, X, Time, r, sigma, table) - + 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, +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] + sigma = table[, 5] } - - # Compute: + + # Compute: I = 0 - xi = (Time*X-I)*exp(-r*Time)/S - (1-exp(-r*Time))/r + 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)) ) + 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, + } 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) +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 + # 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, + # 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) - + # 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] + 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: + + # PDE Parameters - Do not change: T0 = 0 - Tout = 1 + Tout = 1 np = 0 Price.by.S = 0 mf = 12 @@ -822,7 +822,7 @@ 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)) * @@ -830,7 +830,7 @@ 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", @@ -859,11 +859,11 @@ as.double(xbkpt), as.double(work), as.integer(iwork), - PACKAGE = "fAsianOptions" + PACKAGE = "fAsianOptions" ) Price[i] = result[[13]]*S[i] } - + # Return Value: Price } @@ -873,37 +873,37 @@ # LAPLACE INVERSION: -gGemanYor = -function(lambda, S = 100, X = 100, Time = 1, r = 0.05, sigma = 0.30, +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)) - + g = rep(complex(real = 0, imaginary = 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) ) { + 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) + mu = sqrt(2*lambda[i] + nu^2) + q = (sigma^2)*X*Time/(4*S) gamma1 = (mu-nu)/2 - gamma2 = (mu+nu)/2 + gamma2 = (mu+nu)/2 # Convergence Parameters: a = gamma1-2 + 1 b = gamma2+1 + a + 1 @@ -911,22 +911,22 @@ # 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]) + 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 + lam = lambda xlab = "lambda" - ylab = "g" + ylab = "g" } else { - lam = Im(lambda) + lam = Im(lambda) xlab = "Im(lambda)" - ylab = "Re(g)" + ylab = "Re(g)" } lambda.min = 4*r/sigma^2 cat("\nmin lambda:", lambda.min, "\n") @@ -936,17 +936,17 @@ xlab = xlab, ylab = ylab) lines(lam, 0*Re(g)) # Convergence Indexes: - mu = sqrt(2*lambda + nu^2) + 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", + 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)) ) ) + lines(x = lambda.min*c(1,1), y = c( min(c(a,b)), max(c(a,b)) ) ) } - + # Return Value: g } @@ -955,16 +955,16 @@ # ------------------------------------------------------------------------------ -GemanYorAsianOption = -function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, +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 @@ -976,12 +976,12 @@ 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, + zi = complex(real = 0, imaginary = x) + zc = complex(real = alpha, imaginary = 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) ) + Re ( exp(zi*h + g2) / (2*pi) ) } # Call: @@ -992,31 +992,31 @@ if (doprint) { cat("\nDelta:", delta) cat("\nS:", S, "X:", X) - cat("\nTime:", Time, "r:", r, "sigma:", sigma, "\n") + cat("\nTime:", Time, "r:", r, "sigma:", sigma, "\n") } i = 1 - I = integrate(fx, lower = 0, upper = delta, + 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, + 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)) - } - + 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 - } - + Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X + Price = Price - Parity + } + # Return Value: option = list( - price = Price, + price = Price, call = match.call() ) class(option) = "option" option @@ -1027,12 +1027,12 @@ # SPECTRAL EXPANSION: -gLinetzky = -function(x, y, tau, nu, ip = 0) +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 + # 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) @@ -1048,63 +1048,63 @@ for (i in 1:length(x)) { p = x[i] if (p == 0) { - result[i] = 0 + result[i] = 0 } else { z = 1/(2*y) kappa = -(nu+3)/2 - mu = complex(real = 0, imag = p/2) + mu = complex(real = 0, imaginary = 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 + 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 ) + 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 + 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 ) ) + 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) + + 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 + phi = atan(Im(g)/Re(g)) + pi logV4 = r*cos(phi) - argV4 = cos(r*sin(phi)) + argV4 = cos(r*sin(phi)) } # Collect all terms: if ( p < 100) { - result[i] = exp(logV1+logV2+logV3)*V4/(8*pi^2) + result[i] = exp(logV1+logV2+logV3)*V4/(8*pi^2) } else { - result[i] = exp(logV1+logV2+logV3+logV4)*argV4/(8*pi^2) + result[i] = exp(logV1+logV2+logV3+logV4)*argV4/(8*pi^2) } # print(c(p, logV5, logV5a, argV5, argV5a)) } } - - # Return Value: - result + + # 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", +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] @@ -1114,76 +1114,76 @@ 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, + 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 + 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]) + 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 + 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) + 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 } + 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] - } + ( 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) + lines(x, F) } } if (doprint) print(cbind(P, P1)) - + # Derive Call/Put Price: - + # Put Price: - Linetzky = exp(-r*Time) * (S/tau) * P - + 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, + price = Linetzky, call = match.call() ) class(option) = "option" option @@ -1195,14 +1195,14 @@ # 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 +# 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, +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 @@ -1211,50 +1211,50 @@ # 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] + 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 + 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 + 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 + 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, + price = Price, call = match.call() ) class(option) = "option" option @@ -1264,8 +1264,8 @@ # ------------------------------------------------------------------------------ -CurranThompsonAsianOption = -function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, +CurranThompsonAsianOption = +function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz @@ -1273,76 +1273,76 @@ # 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 + # 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) ) + 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, + 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 + 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 + f } - + # Integrate: - value = integrate(f, lower=0, upper=1, g.star=g.star, - alpha=alpha, sigma=sigma, subdivisions=1000, + 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)) ) - + 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 + Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X + CurranThompson = CurranThompson - Parity } - - # Return Value: + + # Return Value: option = list( - price = CurranThompson, + price = CurranThompson, call = match.call() ) class(option) = "option" option } - + # ------------------------------------------------------------------------------ -RogerShiThompsonAsianOption = -function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, +RogerShiThompsonAsianOption = +function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, sigma = 0.30) { # A function implemented by Diethelm Wuertz @@ -1350,165 +1350,165 @@ # 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 + # 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 + # 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) { + gamma.star = function(S, X, r, sigma, lower = -99, upper = 99) { func = function(gamma, S, X, r, sigma) { - f = function(x, gamma, roh, sigma) { + f = function(x, gamma, roh, sigma) { time = x alpha = roh - sigma*sigma/2 - f = exp( 3 * gamma * sigma * time * (1-time/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, + 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 + 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 + 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 + f } - + # Integrate: - value = integrate(f, lower=0, upper=1, g.star=g.star, roh=r, + 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)) ) - + 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 + Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X RogerShiThompson = RogerShiThompson - Parity } - - # Return Value: + + # Return Value: option = list( - price = RogerShiThompson, + price = RogerShiThompson, call = match.call() ) class(option) = "option" option } - + # ------------------------------------------------------------------------------ -ThompsonAsianOption = -function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, +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. - + # 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 + # 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, + 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 = (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 } + 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 } + 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)) } - + 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, + 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 - + 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 + Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X Thompson = Thompson - Parity } - - # Return Value: + + # Return Value: option = list( - price = Thompson, + price = Thompson, call = match.call() ) class(option) = "option" option @@ -1518,27 +1518,27 @@ # ------------------------------------------------------------------------------ -TolmatzAsianOption = -function(TypeFlag = c("c", "p"), S = 100, X = 100, Time = 1, r = 0.09, +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 - + # 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 + Parity = (1/(r*Time))*(1-exp(-r*Time))*S - exp(-r*Time)*X Tolmatz = Tolmatz - Parity } - - # Return Value: + + # Return Value: option = list( - price = Tolmatz, + price = Tolmatz, call = match.call() ) class(option) = "option" option @@ -1549,8 +1549,8 @@ # SYMMETRY AND EQUIVALENCE RELATIONS: -CallPutParityAsianOption = -function(TypeFlag = "p", Price = 8.828759, S = 100, X = 100, Time = 1, +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 @@ -1561,17 +1561,17 @@ 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) + Parity = S*(1-exp(-r*Time))/(r*Time) - X*exp(-r*Time) result = Price - Parity } # Return Value: @@ -1582,53 +1582,53 @@ # ------------------------------------------------------------------------------ -WithDividendsAsianOption = -function(TypeFlag = "c", Dividends = 0.45, S = 100, X = 100, Time = 1, +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, + 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, + price = Price, call = match.call() ) class(option) = "option" - option -} + option +} ################################################################################ # TABULATED RESULTS: -FuMadanWangTable = +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, @@ -1646,21 +1646,21 @@ 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) + 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) + 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, 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, + S, X, r, sigma, CallEulerFMW, CPUEulerFMW, CallPostWidderFMW, CPUPostWidderFMW, CallZhang)) } @@ -1669,24 +1669,24 @@ # ------------------------------------------------------------------------------ -FusaiTaglianiTable = +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, @@ -1700,8 +1700,8 @@ 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) - + CallFTLN = c(CallLNa, CallLNb, CallLNc) + # Gram-Charlier Call Prices: CallFTa = c( 11.95127, 7.40754, 3.64091, 1.31097, 0.33156, @@ -1716,7 +1716,7 @@ 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) } @@ -1725,31 +1725,31 @@ # ------------------------------------------------------------------------------ -GemanTable = +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)) } @@ -1758,41 +1758,41 @@ # ------------------------------------------------------------------------------ -LinetzkyTable = +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 - + 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.0559860415, 0.2183875466, 0.1722687410, 0.1931737903, 0.2464156905, 0.3062203648, 0.3500952190) CallSLT = c( - 0.055986, 0.218388, 0.172269, 0.193174, + 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.05602, 0.2185, 0.1725, 0.1933, 0.2465, 0.3064, 0.3503) CallTLB = c( - 0.055985, 0.218366, 0.172226, 0.193060, + 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.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) } @@ -1801,25 +1801,25 @@ # ------------------------------------------------------------------------------ -ZhangTable = +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, @@ -1828,7 +1828,7 @@ 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)) } @@ -1837,25 +1837,25 @@ # ------------------------------------------------------------------------------ -ZhangShortTable = +ZhangShortTable = function() { # A function implemented by Diethelm Wuertz - + # Description: # Display Zhangs's short-tenor (ST) results from Table 7 - + # Source: - # J.E. Zhang [2002], + # 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, @@ -1865,7 +1865,7 @@ 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)) } @@ -1874,25 +1874,25 @@ # ------------------------------------------------------------------------------ -ZhangLongTable = +ZhangLongTable = function() { # A function implemented by Diethelm Wuertz - + # Description: # Display Zhangs's long-tenor (LT) results from Table 7 - + # Source: - # J.E. Zhang [2002], + # 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, @@ -1902,7 +1902,7 @@ 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)) } diff -Nru fasianoptions-2100.76/R/GammaFunctions.R fasianoptions-3000.78/R/GammaFunctions.R --- fasianoptions-2100.76/R/GammaFunctions.R 2009-09-02 20:27:40.000000000 +0000 +++ fasianoptions-3000.78/R/GammaFunctions.R 2013-01-11 11:21:37.000000000 +0000 @@ -6,31 +6,31 @@ # # 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 +# 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, +# 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) -# this R-port: +# Copyrights (C) +# this R-port: # by Diethelm Wuertz # for the code accessed (or partly included) from other R-ports: # R: see R's copyright and license file # for Haug's Option Pricing Formulas: -# Formulas are implemented along the book and the Excel spreadsheets of +# Formulas are implemented along the book and the Excel spreadsheets of # E.G. Haug, "The Complete Guide to Option Pricing"; documentation # is partly taken from www.derivicom.com which implements -# a C Library based on Haug. For non-academic and commercial use -# we recommend the professional software from "www.derivicom.com". +# a C Library based on Haug. For non-academic and commercial use +# we recommend the professional software from "www.derivicom.com". ################################################################################ # FUNCTION: DESCRIPTION: # erf Error function -# [gamma] Gamma function +# [gamma] Gamma function # [lgamma] LogGamma function, returns log(gamma) # [digamma] First derivative of of LogGamma, dlog(gamma(x))/dx # [trigamma] Second derivative of of LogGamma, dlog(gamma(x))/dx @@ -49,22 +49,22 @@ ################################################################################ -erf = +erf = function(x) { # A function implemented by Diethelm Wuertz # Description: # Computes the Error function for real argument "x" - + # Arguments: # x - a real numeric value or vector. - + # FUNCTION: - + # Result # DW 2005-05-04 ans = 2 * pnorm(sqrt(2) * x) - 1 - + # Return Value: ans } @@ -73,45 +73,45 @@ # ------------------------------------------------------------------------------ -cgamma = -function(x, log = FALSE) +cgamma = +function(x, log = FALSE) { # A function implemented by Diethelm Wuertz # Description: - # Computes the Gamma Function for complex argument "x" - + # Computes the Gamma Function for complex argument "x" + # Arguments: # z - a complex or real vector # log - if TRUE the logarithm of the gamma is calculated - # otherwise if FALSE, the gamma function itself + # otherwise if FALSE, the gamma function itself # will be calculated. - + # Source: # For the Fortran Routine: # http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html - + # FUNCTION: - + # Test for complex arguments: - if (!is.complex(x)) x = complex(real = x, imag = 0*x) - + if (!is.complex(x)) x = complex(real = x, imaginary = 0*x) + # Calculate Gamma: KF = 1 if (log) { KF = KF - 1 } result = rep(NA, times = length(x)) - for ( i in 1:length(x) ) { + for ( i in 1:length(x) ) { value = .Fortran("cgama", as.double(Re(x[i])), as.double(Im(x[i])), as.integer(KF), as.double(0), - as.double(0), + as.double(0), PACKAGE = "fAsianOptions") - result[i] = complex(real = value[[4]], imag = value[[5]]) + result[i] = complex(real = value[[4]], imaginary = value[[5]]) } - + # Return Value: result } @@ -120,46 +120,46 @@ # ------------------------------------------------------------------------------ -Psi = +Psi = function(x) { # A function implemented by Diethelm Wuertz # Description: # Computes the Psi or Digamma function for complex or real argument - + # Arguments: # z - a complex numeric value or vector. - + # Details: # [AS} formula 6.3.1 # $ \Psi(x) = d ln \Gamma(z) / dz = \Gamma prime (z) / \Gamma(z) $ - + # Arguments: # x - complex or real vector - + # Source: # For the Fortran Routine: # http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html - + # FUNCTION: - + # Psi: result = rep(NA, times = length(x)) if (!is.complex(x)) { # Use R's digamma() function: - result = digamma(x) + result = digamma(x) } else { - for ( i in 1:length(Re(x)) ) { + for ( i in 1:length(Re(x)) ) { value = .Fortran("cpsi", as.double(Re(x[i])), as.double(Im(x[i])), as.double(0), as.double(0), PACKAGE = "fAsianOptions") - result[i] = complex(real = value[[3]], imag = value[[4]]) + result[i] = complex(real = value[[3]], imaginary = value[[4]]) } } - + # Return Value: result } @@ -168,41 +168,41 @@ # ------------------------------------------------------------------------------ -igamma = +igamma = function(x, a) { # A function implemented by Diethelm Wuertz # Description: - # Computes the Incomplete Gamma Function P(a, x) with - # Re(a) > 0 for complex or real argument "x" and for - # complex or real index "z" - + # Computes the Incomplete Gamma Function P(a, x) with + # Re(a) > 0 for complex or real argument "x" and for + # complex or real index "z" + # Arguments: # z - a complex or real vector # a - a complex or real numeric value - + # Details: - # [AS] formula 6.5.1 + # [AS] formula 6.5.1 # $ frac{1}{Gamma(a)} * \int_0^x e^{-t} t^{a-1} dt $ - + # FUNCTION: - + # igamma: if (!is.complex(x) && !is.complex(a)) { # Use R's pgamma() function: # if (a < 0) Not suppported ... - result = pgamma(x, a) + result = pgamma(x, a) } else { # Why not derive the result from KummersM ? log = FALSE if (log) { # Not yet supported: - result = kummerM(a, a + 1, -x, lnchf = 1) + a*log(x) - log(a) + result = kummerM(a, a + 1, -x, lnchf = 1) + a*log(x) - log(a) } else { - result = kummerM(a, a + 1, -x, lnchf = 0) * x^a / a - } + result = kummerM(a, a + 1, -x, lnchf = 0) * x^a / a + } } - + # Return Value: result } @@ -211,47 +211,48 @@ # ------------------------------------------------------------------------------ -Pochhammer = +Pochhammer = function(x, n) { # A function implemented by Diethelm Wuertz - + # Description: # Computes Pochhammer's Symbol - + # Arguments: # x - a complex numeric value or vector. - # n - an integer n >=0. An notation used in the theory of special - # functions for the rising factorial, also known as the rising - # factorial power (Graham et al. 1994). - + # n - an integer n >=0. An notation used in the theory of special + # functions for the rising factorial, also known as the rising + # factorial power (Graham et al. 1994). + # Details: # as defined in [AS] by formula 6.1.22 - + # FUNCTION: - + # Note: # $ (z)_0 = 1 $ # $ (z)_n = z(z+1)(z+2) \dots (z+n-1) = frac{\Gamma(z+n)}{Gamma(z)} $ - + # In case of wrong argument Type: Pochhammer = NA - + # For Complex Arguments: if (is.complex(x)) { - Pochhammer = cgamma(x + n)/cgamma(x) + Pochhammer = cgamma(x + n)/cgamma(x) } - + # For Real Arguments: # DW: 2006-05-10 is.real(z) replaced by is.real(x) - if (is.real(x)) { - Pochhammer = gamma(x + n)/gamma(x) + # YC: is.real is deprecated -> replaced by is.double + if (is.double(x)) { + Pochhammer = gamma(x + n)/gamma(x) } - + # Return Value: Pochhammer } - + ################################################################################ # SPlus Addon for beta() and lbeta() # quick and dirty implementation ... @@ -267,33 +268,32 @@ beta = function(a, b) { # A function implemented by Diethelm Wuertz - + # Description: # Computes the beta function - + # Result: - ans = gamma(a) * gamma(b) / gamma(a+b) - + ans = gamma(a) * gamma(b) / gamma(a+b) + # Return Value: - ans + ans } - - - lbeta = + + + lbeta = function(a, b) { # A function implemented by Diethelm Wuertz - + # Description: # Computes the beta function - + # Result: ans = lgamma(a) + lgamma(b) - lgamma(a+b) - + # Return Value: - ans + ans } } ################################################################################ - diff -Nru fasianoptions-2100.76/R/HypergeometricFunctions.R fasianoptions-3000.78/R/HypergeometricFunctions.R --- fasianoptions-2100.76/R/HypergeometricFunctions.R 2009-09-02 20:27:40.000000000 +0000 +++ fasianoptions-3000.78/R/HypergeometricFunctions.R 2012-03-19 10:09:48.000000000 +0000 @@ -15,7 +15,7 @@ # MA 02111-1307 USA # Copyrights (C) -# for this R-port: +# for this R-port: # 1999 - 2004, Diethelm Wuertz, GPL # Diethelm Wuertz # info@rmetrics.org @@ -28,7 +28,7 @@ ################################################################################ -# FUNCTION: KUMMER DESCRIPTION: +# FUNCTION: KUMMER DESCRIPTION: # kummerM Computes Confluent Hypergeometric Function of the 1st Kind # kummerU Computes Confluent Hypergeometric Function of the 2nd Kind # FUNCTION: WHITTAKER DESCRIPTION: @@ -43,29 +43,29 @@ # KUMMER: -kummerM = -function(x, a, b, lnchf = 0, ip = 0) +kummerM = +function(x, a, b, lnchf = 0, ip = 0) { # A function implemented by Diethelm Wuertz # Description: - # Calculate the Confluent Hypergeometric Function of the First + # Calculate the Confluent Hypergeometric Function of the First # Kind for complex argument "x" and complex indexes "a" and "b" - + # Arguments: # x - complex function argument - # a, b - complex indexes - # lnchf - - # ip - - + # a, b - complex indexes + # lnchf - + # ip - + # FUNCTION: - + # You can also input real arguments: - if (!is.complex(x)) x = complex(real = x, imag = 0*x) - if (!is.complex(a)) a = complex(real = a, imag = 0) - if (!is.complex(b)) b = complex(real = b, imag = 0) - + if (!is.complex(x)) x = complex(real = x, imaginary = 0*x) + if (!is.complex(a)) a = complex(real = a, imaginary = 0) + if (!is.complex(b)) b = complex(real = b, imaginary = 0) + # Calculate KummerM: - chm = rep(complex(real = 0, imag = 0), length = length(x)) + chm = rep(complex(real = 0, imaginary = 0), length = length(x)) value = .Fortran("chfm", as.double(Re(x)), as.double(Im(x)), @@ -79,8 +79,8 @@ as.integer(lnchf), as.integer(ip), PACKAGE = "fAsianOptions") - result = complex(real = value[[7]], imag = value[[8]]) - + result = complex(real = value[[7]], imaginary = value[[8]]) + # Return Value: result } @@ -89,38 +89,38 @@ # ------------------------------------------------------------------------------ -kummerU = -function(x, a, b, ip = 0) +kummerU = +function(x, a, b, ip = 0) { # A function implemented by Diethelm Wuertz # Description: - # Calculate the Confluent Hypergeometric Function of the Second + # Calculate the Confluent Hypergeometric Function of the Second # Kind for complex argument "x" and complex indexes "a" and "b" - + # Arguments: - + # FUNCTION: - + # Todo ... lnchf = 0 - + # Test for complex arguments: - if (!is.complex(x)) x = complex(real = x, imag = 0*x) - if (!is.complex(a)) a = complex(real = a, imag = 0) - if (!is.complex(b)) b = complex(real = b, imag = 0) - + if (!is.complex(x)) x = complex(real = x, imaginary = 0*x) + if (!is.complex(a)) a = complex(real = a, imaginary = 0) + if (!is.complex(b)) b = complex(real = b, imaginary = 0) + # Calculate KummerU: # From KummerM: # Uses the formula ... - # pi/sin(pi*b) [ M(a,b,z) / (Gamma(1+a-b)*Gamma(b)) - + # pi/sin(pi*b) [ M(a,b,z) / (Gamma(1+a-b)*Gamma(b)) - # x^(1-b) * M(1+a-b,2-b,z) / (Gamma(a)*Gamma(2-b)) ] ans = ( pi/sin(pi*b) ) * ( kummerM(x, a = a, b = b, lnchf = lnchf, ip=ip) / - ( cgamma(1+a-b)*cgamma(b) ) - (x^(1-b)) * - kummerM(x, a = (1+a-b), b=2-b, lnchf = lnchf, ip = ip) / + ( cgamma(1+a-b)*cgamma(b) ) - (x^(1-b)) * + kummerM(x, a = (1+a-b), b=2-b, lnchf = lnchf, ip = ip) / ( cgamma(a)*cgamma(2-b) ) ) - # Return Value: + # Return Value: ans } @@ -129,25 +129,25 @@ # WHITTAKER: -whittakerM = -function(x, kappa, mu, ip = 0) +whittakerM = +function(x, kappa, mu, ip = 0) { # A function implemented by Diethelm Wuertz # Description: # Computes Whittaker's M Function - + # Arguments: - + # FUNCTION: - + # Test for complex arguments: - if (!is.complex(x)) x = complex(real = x, imag = 0*x) - if (!is.complex(kappa)) kappa = complex(real = kappa, imag = 0) - if (!is.complex(mu)) mu = complex(real = mu, imag = 0) - + if (!is.complex(x)) x = complex(real = x, imaginary = 0*x) + if (!is.complex(kappa)) kappa = complex(real = kappa, imaginary = 0) + if (!is.complex(mu)) mu = complex(real = mu, imaginary = 0) + # Calculate: ans = exp(-x/2) * x^(1/2+mu) * kummerM(x, 1/2+mu-kappa, 1+2*mu, ip = ip) - + # Return Value: ans } @@ -156,25 +156,25 @@ # ------------------------------------------------------------------------------ -whittakerW = -function(x, kappa, mu, ip = 0) +whittakerW = +function(x, kappa, mu, ip = 0) { # A function implemented by Diethelm Wuertz # Description: # Computes Whittaker's M Function - + # Arguments: - + # FUNCTION: - + # Test for complex arguments: - if (!is.complex(x)) x = complex(real = x, imag = 0*x) - if (!is.complex(kappa)) kappa = complex(real = kappa, imag = 0) - if (!is.complex(mu)) mu = complex(real = mu, imag = 0) - + if (!is.complex(x)) x = complex(real = x, imaginary = 0*x) + if (!is.complex(kappa)) kappa = complex(real = kappa, imaginary = 0) + if (!is.complex(mu)) mu = complex(real = mu, imaginary = 0) + # Calculate: ans = exp(-x/2) * x^(1/2+mu) * kummerU(x, 1/2+mu-kappa, 1+2*mu, ip = ip) - + # Return Value: ans } @@ -184,27 +184,27 @@ # HERMITE POLYNOMIAL: -hermiteH = +hermiteH = function(x, n, ip = 0) { # A function implemented by Diethelm Wuertz - + # Description: - # Computes the Hermite Polynomial - + # Computes the Hermite Polynomial + # Arguments: # n - the index of the Hermite polynomial. - + # FUNCTION: - + # Check stopifnot(n - round(n, 0) == 0) - + # Result: S = sign(x) + (1-sign(abs(x))) ans = (S*2)^n * Re ( kummerU(x^2, -n/2, 1/2, ip = ip) ) - + # Return Value: - ans + ans } diff -Nru fasianoptions-2100.76/R/zzz.R fasianoptions-3000.78/R/zzz.R --- fasianoptions-2100.76/R/zzz.R 2009-09-02 20:27:40.000000000 +0000 +++ fasianoptions-3000.78/R/zzz.R 2013-01-11 12:24:51.000000000 +0000 @@ -30,26 +30,9 @@ ################################################################################ -.First.lib = -function(lib, pkg) -{ - # Startup Mesage and Desription: - MSG <- if(getRversion() >= "2.5") packageStartupMessage else message - dsc <- packageDescription(pkg) - if(interactive() || getOption("verbose")) { - # not in test scripts - MSG(sprintf("Rmetrics Package %s (%s) loaded.", pkg, dsc$Version)) - } - - # Load dll: - library.dynam("fAsianOptions", pkg, lib) -} - - if(!exists("Sys.setenv", mode = "function")) # pre R-2.5.0, use "old form" Sys.setenv <- Sys.putenv ################################################################################ - diff -Nru fasianoptions-2100.76/debian/changelog fasianoptions-3000.78/debian/changelog --- fasianoptions-2100.76/debian/changelog 2013-05-05 16:25:25.000000000 +0000 +++ fasianoptions-3000.78/debian/changelog 2013-05-05 16:25:25.000000000 +0000 @@ -1,3 +1,36 @@ +fasianoptions (3000.78-2precise0) precise; urgency=low + + * Compilation for Ubuntu 12.04.2 LTS + + -- Michael Rutter Sun, 05 May 2013 16:16:10 +0000 + +fasianoptions (3000.78-2) unstable; urgency=low + + * debian/control: Set Build-Depends: to current R version + + * (Re-)building with R 3.0.0 + + -- Dirk Eddelbuettel Wed, 03 Apr 2013 05:45:07 -0500 + +fasianoptions (3000.78-1) unstable; urgency=low + + * New upstream release + + * debian/control: Set Build-Depends: to current R version + * debian/control: Set Standards-Version: to current version + + -- Dirk Eddelbuettel Fri, 08 Feb 2013 13:37:23 -0600 + +fasianoptions (2160.77-1) unstable; urgency=low + + * New upstream release (Closes: #665235) + + * debian/control: Set Build-Depends: to current R version + * debian/control: Change Depends to ${R:Depends} + * debian/control: Set Standards-Version: to current version + + -- Dirk Eddelbuettel Thu, 22 Mar 2012 12:00:43 -0500 + fasianoptions (2100.76-3) unstable; urgency=low * Rebuilt for R 2.14.0 so that a default NAMESPACE file is created diff -Nru fasianoptions-2100.76/debian/control fasianoptions-3000.78/debian/control --- fasianoptions-2100.76/debian/control 2013-05-05 16:25:25.000000000 +0000 +++ fasianoptions-3000.78/debian/control 2013-05-05 16:25:25.000000000 +0000 @@ -2,13 +2,13 @@ Section: gnu-r Priority: optional Maintainer: Dirk Eddelbuettel -Build-Depends: debhelper (>= 7.0.0), r-base-dev (>> 2.13.2), cdbs, r-cran-fbasics (>= 290.76), r-cran-foptions (>= 270.74), xvfb, xauth, xfonts-base -Standards-Version: 3.9.2 +Build-Depends: debhelper (>= 7.0.0), r-base-dev (>= 3.0.0~), cdbs, r-cran-fbasics (>= 2160.81-2), r-cran-foptions (>= 270.74), xvfb, xauth, xfonts-base +Standards-Version: 3.9.4 Homepage: http://www.Rmetrics.org Package: r-cran-fasianoptions Architecture: any -Depends: ${shlibs:Depends}, r-base-core (>> 2.13.2), r-cran-foptions (>= 270.74), r-cran-fbasics (>= 290.76) +Depends: ${shlibs:Depends}, ${R:Depends}, r-cran-foptions (>= 270.74), r-cran-fbasics (>= 2160.81-2) Suggests: r-cran-runit Description: GNU R package for financial engineering -- fAsianOptions This package of functions for financial engineering and computational Binary files /tmp/53AimPR2_W/fasianoptions-2100.76/inst/DocCopying.pdf and /tmp/ByVZSaEYbf/fasianoptions-3000.78/inst/DocCopying.pdf differ diff -Nru fasianoptions-2100.76/src/EBMAsianOptions.f fasianoptions-3000.78/src/EBMAsianOptions.f --- fasianoptions-2100.76/src/EBMAsianOptions.f 2009-09-02 20:27:40.000000000 +0000 +++ fasianoptions-3000.78/src/EBMAsianOptions.f 2013-02-08 15:12:54.000000000 +0000 @@ -1,138 +1,146 @@ +C added R print function to ensure that output goes to console output + +C ------------------------------------------------------------------------------ + C ALGORITHM 540R (REMARK ON ALG.540), COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 343-344. C C USED FOR: C PROGRAM RUNPDE -C FOR TESTING AND DEBUGGING UNDER FORTRAN -C CALL PDETEST() +C FOR TESTING AND DEBUGGING UNDER FORTRAN +C CALL PDETEST() C END C MODEL 1: VECER's PDE C MODEL 2: ZHANG's PDE +C Small modification to avoid gfortran warning: 'Rank mismatch in +C argument 'x' at (1) (rank-1 and scalar)' by Yohan Chalabi on +C March 2012 + C ------------------------------------------------------------------------------ +c$$$ +c$$$ SUBROUTINE PDETEST() +c$$$CC NOT USED BY R +c$$$CC FOR TESTING AND DEBUGGING UNDER FORTRAN +c$$$ IMPLICIT REAL*8 (A-H, O-Z) +c$$$ PARAMETER(MNP=10) +c$$$ DIMENSION PRICE(MNP+1), XBYS(MNP+1) +c$$$C +c$$$ PARAMETER (MMF=12, MMX=1000) +c$$$ PARAMETER (MNPDE=1, MKORD=4, MNINT=MMX, MNCC=2, MMAXDER=5) +c$$$ +c$$$C WORKING ARRAYS: +c$$$ DIMENSION WORK +c$$$ * (MKORD+MNPDE*(4+9*MNPDE)+(MKORD+(MNINT-1)*(MKORD-MNCC))* +c$$$ * (3*MKORD+2+MNPDE*(3*(MKORD-1)*MNPDE+MMAXDER+4))) +c$$$ DIMENSION IWORK((MNPDE+1)*(MKORD+(MNINT-1)*(MKORD-MNCC))) +c$$$ DIMENSION XBKPT(MNINT+1) +c$$$ +c$$$C PDE PARAMETERS: +c$$$ NP = MNP +c$$$ MF = MMF +c$$$ MX = MMX +c$$$ NPDE = MNPDE +c$$$ KORD = MKORD +c$$$ NINT = MNINT +c$$$ NCC = MNCC +c$$$ MAXDER = MMAXDER +c$$$ +c$$$C OPTION SETTINGS: +c$$$ SIGMA = 0.30D0 +c$$$ TIME = 1.00D0 +c$$$ RR = 0.09D0 +c$$$ XS = 100.00D0 +c$$$ XSMIN = 90.00D0 +c$$$ XSMAX = 110.00D0 +c$$$ SS = 100.00D0 +c$$$ DELTA = (XSMAX-XSMIN)/NP +c$$$ DO I = 1, NP+1 +c$$$ XBYS(I) = (XSMIN +(I-1)*DELTA)/XS +c$$$ ENDDO +c$$$ +c$$$C SET TIME POINTS: +c$$$C T0 = INITIAL VALUE OF T, THE INDEPENDENT VARIABLE +c$$$C TOUT = VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT +c$$$C DT = INITIAL STEP SIZE IN T +c$$$C EPS = RELATIVE TIME ERROR BOUND +c$$$ T0 = 0.0D0 +c$$$ TOUT = 1.0D0 +c$$$ EPS = 1.0D-08 +c$$$ DT = 1.0D-10 +c$$$ +c$$$C FURTHER PARAMETERS: +c$$$C NINT=1000 - NUMBER OF SUBINTERVALS (XLEFT,XRIGHT) IS TO BE DIVIDED +c$$$C KORD=4 - ORDER OF THE PIECEWISE POLYNOMIAL SPACE TO BE USED +c$$$C NCC=2 - NUMBER OF CONTINUITY CONDITIONS TO BE IMPOSED +c$$$C MF=12 - METHOD FLAG +c$$$C ADAMS METHODS - GENERALIZATIONS OF CRANK-NICOLSON AND +c$$$C CHORD METHOD WITH FINITE DIFFERENCES JACOBIAN +c$$$C INDEX - INTEGER USED ON INPUT TO INDICATE TYPE OF CALL +c$$$C WORK - WORKING ARRAY +c$$$C IWORK - SIZE OF WORKING ARRAY +c$$$ +c$$$C ASIAN CALL (1) AND PUT(2) VALUE: +c$$$ Z = -1 +c$$$ DO IP = 1, 2 +c$$$ Z = -Z +c$$$C PDE PARAMETERS: +c$$$ MODSEL = 1 +c$$$ SIGMAT = SIGMA * DSQRT(TIME) +c$$$ RRT = RR*TIME +c$$$ XM = 5.0D0 * SIGMAT +c$$$ WRITE (*,*) +c$$$ WRITE (*,*) " PDE - ASIAN OPTION SETTINGS" +c$$$ WRITE (*,*) " MF KORD NCC : ", MF, KORD, NCC +c$$$ WRITE (*,*) " SIGMA*TIME : ", SIGMAT +c$$$ WRITE (*,*) " R*TIME : ", RRT +c$$$ WRITE (*,*) " XM : ", XM +c$$$ WRITE (*,*) " (XMIN,XMAX)/S : ", XSMIN/SS, XSMAX/SS +c$$$ CALL ASIANVAL( +c$$$ & Z, SS, XS, XSMIN, XSMAX, TIME, RR, SIGMA, +c$$$ & T0, TOUT, EPS, DT, PRICE, NP, MODSEL, +c$$$ & MF, NPDE, KORD, MX, NCC, MAXDER, +c$$$ & XBYS, XBKPT, WORK, IWORK) +c$$$C OUTPUT U - NUMERICAL SOLUTION: +c$$$ WRITE (*,*) " XLEFT XRIGHT : ", XBKPT(1), XBKPT(NINT+1) +c$$$ WRITE (*,*) " EPS DT MX : ", EPS, DT, MX +c$$$ WRITE (*,*) " ERROR CODE: : ", INDEX +c$$$ WRITE(*,*) +c$$$ WRITE(*,*) " U - NUMERICAL SOLUTION FOR DIFF STRIKES:" +c$$$ WRITE(*,*) " X XI PRICE " +c$$$ DO I = 1, NP+1 +c$$$ XI = XBYS(I)*EXP(-RRT) - (1.0-EXP(-RRT))/RRT +c$$$ WRITE(*,9) XS*XBYS(I), SS*XI, SS*PRICE(I), SS*(PRICE(I)-XI) +c$$$ ENDDO +c$$$ ENDDO +c$$$ 9 FORMAT(F10.3, 4F15.7) +c$$$ +c$$$ RETURN +c$$$ END - SUBROUTINE PDETEST() -CC NOT USED BY R -CC FOR TESTING AND DEBUGGING UNDER FORTRAN - IMPLICIT REAL*8 (A-H, O-Z) - PARAMETER(MNP=10) - DIMENSION PRICE(MNP+1), XBYS(MNP+1) -C - PARAMETER (MMF=12, MMX=1000) - PARAMETER (MNPDE=1, MKORD=4, MNINT=MMX, MNCC=2, MMAXDER=5) -C WORKING ARRAYS: - DIMENSION WORK - * (MKORD+MNPDE*(4+9*MNPDE)+(MKORD+(MNINT-1)*(MKORD-MNCC))* - * (3*MKORD+2+MNPDE*(3*(MKORD-1)*MNPDE+MMAXDER+4))) - DIMENSION IWORK((MNPDE+1)*(MKORD+(MNINT-1)*(MKORD-MNCC))) - DIMENSION XBKPT(MNINT+1) - -C PDE PARAMETERS: - NP = MNP - MF = MMF - MX = MMX - NPDE = MNPDE - KORD = MKORD - NINT = MNINT - NCC = MNCC - MAXDER = MMAXDER - -C OPTION SETTINGS: - SIGMA = 0.30D0 - TIME = 1.00D0 - RR = 0.09D0 - XS = 100.00D0 - XSMIN = 90.00D0 - XSMAX = 110.00D0 - SS = 100.00D0 - DELTA = (XSMAX-XSMIN)/NP - DO I = 1, NP+1 - XBYS(I) = (XSMIN +(I-1)*DELTA)/XS - ENDDO - -C SET TIME POINTS: -C T0 = INITIAL VALUE OF T, THE INDEPENDENT VARIABLE -C TOUT = VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT -C DT = INITIAL STEP SIZE IN T -C EPS = RELATIVE TIME ERROR BOUND - T0 = 0.0D0 - TOUT = 1.0D0 - EPS = 1.0D-08 - DT = 1.0D-10 - -C FURTHER PARAMETERS: -C NINT=1000 - NUMBER OF SUBINTERVALS (XLEFT,XRIGHT) IS TO BE DIVIDED -C KORD=4 - ORDER OF THE PIECEWISE POLYNOMIAL SPACE TO BE USED -C NCC=2 - NUMBER OF CONTINUITY CONDITIONS TO BE IMPOSED -C MF=12 - METHOD FLAG -C ADAMS METHODS - GENERALIZATIONS OF CRANK-NICOLSON AND -C CHORD METHOD WITH FINITE DIFFERENCES JACOBIAN -C INDEX - INTEGER USED ON INPUT TO INDICATE TYPE OF CALL -C WORK - WORKING ARRAY -C IWORK - SIZE OF WORKING ARRAY - -C ASIAN CALL (1) AND PUT(2) VALUE: - Z = -1 - DO IP = 1, 2 - Z = -Z -C PDE PARAMETERS: - MODSEL = 1 - SIGMAT = SIGMA * DSQRT(TIME) - RRT = RR*TIME - XM = 5.0D0 * SIGMAT - WRITE (*,*) - WRITE (*,*) " PDE - ASIAN OPTION SETTINGS" - WRITE (*,*) " MF KORD NCC : ", MF, KORD, NCC - WRITE (*,*) " SIGMA*TIME : ", SIGMAT - WRITE (*,*) " R*TIME : ", RRT - WRITE (*,*) " XM : ", XM - WRITE (*,*) " (XMIN,XMAX)/S : ", XSMIN/SS, XSMAX/SS - CALL ASIANVAL( - & Z, SS, XS, XSMIN, XSMAX, TIME, RR, SIGMA, - & T0, TOUT, EPS, DT, PRICE, NP, MODSEL, - & MF, NPDE, KORD, MX, NCC, MAXDER, - & XBYS, XBKPT, WORK, IWORK) -C OUTPUT U - NUMERICAL SOLUTION: - WRITE (*,*) " XLEFT XRIGHT : ", XBKPT(1), XBKPT(NINT+1) - WRITE (*,*) " EPS DT MX : ", EPS, DT, MX - WRITE (*,*) " ERROR CODE: : ", INDEX - WRITE(*,*) - WRITE(*,*) " U - NUMERICAL SOLUTION FOR DIFF STRIKES:" - WRITE(*,*) " X XI PRICE " - DO I = 1, NP+1 - XI = XBYS(I)*EXP(-RRT) - (1.0-EXP(-RRT))/RRT - WRITE(*,9) XS*XBYS(I), SS*XI, SS*PRICE(I), SS*(PRICE(I)-XI) - ENDDO - ENDDO - 9 FORMAT(F10.3, 4F15.7) - - RETURN - END - - C ------------------------------------------------------------------------------ - + SUBROUTINE ASIANVAL( - & ZZ, SS1, XS1, XSMIN, XSMAX, TIME1, RR1, SIGMA1, + & ZZ, SS1, XS1, XSMIN, XSMAX, TIME1, RR1, SIGMA1, & T0, TOUT, EPS, DT, PRICEBYS, NP, MODSEL, - & MF1, NPDE1, KORD1, MX1, NCC1, MAXDER1, + & MF1, NPDE1, KORD1, MX1, NCC1, MAXDER1, & XBYS, XBKPT, WORK, IWORK) - + IMPLICIT REAL*8 (A-H, O-Z) - PARAMETER(MKORD=4, MDERV=0) + PARAMETER(MKORD=4, MDERV=0) DIMENSION WORK * (KORD1+NPDE1*(4+9*NPDE1)+(KORD1+(MX1-1)*(KORD1-NCC1))* * (3*KORD1+2+NPDE1*(3*(KORD1-1)*NPDE1+MAXDER1+4))) DIMENSION IWORK((NPDE1+1)*(KORD1+(MX1-1)*(KORD1-NCC1))) - DIMENSION XBKPT(MX1+1) + DIMENSION XBKPT(MX1+1) DIMENSION USOL(1,1,MDERV+1), SCRTCH(MKORD*(MDERV+1)) DIMENSION XBYS(NP), PRICEBYS(NP) - + COMMON /SIZES/ NINT,KORD,NCC,NPDE,NCPTS,NEQN,IQUAD COMMON /GEAR0/ HUSED, NQUSED, NS, NF, NJ COMMON /GEAR1/ T,DTC,DTMN,DTMX,EPSC,UROUND,N,MFC,KFLAG,JSTART @@ -143,6 +151,8 @@ COMMON /ASIAN1/ SIGMAT, RRT, XM, Z, MODEL COMMON /ASIAN2/ SIGMA, TIME, RR, XS, SS, ETA, XL, XR + DIMENSION XI(1) + C FOR COMMON BLOCKS: SIGMA = SIGMA1 TIME = TIME1 @@ -150,7 +160,7 @@ XS = XS1 SS = SS1 -C FOR COMMON BLOCKS: +C FOR COMMON BLOCKS: MF = MF1 NPDE = NPDE1 KORD = KORD1 @@ -159,78 +169,78 @@ MAXDER = MAXDER1 NINT = MX1 MODEL = MODSEL - PI = 4.0D0 * DATAN(1.0D0) + PI = 4.0D0 * DATAN(1.0D0) -C CALCULATE FOR BOTH, FOR A CALL Z=+1 OR FOR A PUT Z=-1: +C CALCULATE FOR BOTH, FOR A CALL Z=+1 OR FOR A PUT Z=-1: Z = ZZ - -C WORKSPACE SETTINGS: + +C WORKSPACE SETTINGS: IWORK(1) = KORD+NPDE*(4+9*NPDE)+(KORD+(MX-1)* * (KORD-NCC))*(3*KORD+2+NPDE*(3*(KORD-1)*NPDE+MAXDER+4)) - IWORK(2) = (NPDE+1)*(KORD+(NINT-1)*(KORD-NCC)) + IWORK(2) = (NPDE+1)*(KORD+(NINT-1)*(KORD-NCC)) DO I = 1, IWORK(1) WORK(I)=0.0 ENDDO - -C OPTION SETTINGS: + +C OPTION SETTINGS: SIGMAT = SIGMA * DSQRT(TIME) - RRT = RR*TIME - XM = 5.0D0 * SIGMAT + RRT = RR*TIME + XM = 5.0D0 * SIGMAT XL = -XM XR = +XM ETA = (SIGMA**2)*(TIME**3)/6.0D0 - -C SET SPACE POINTS: + +C SET SPACE POINTS: NX = MX DX = 2.0D0 * XM / NX DO I = 1, NX + 1 XBKPT(I) = -XM + (I-1)*DX ENDDO - + C SOLVE PDE: INDEX = 1 - CALL PDECOL(T0, TOUT, DT, XBKPT, EPS, + CALL PDECOL(T0, TOUT, DT, XBKPT, EPS, & NX, KORD, NCC, NPDE, MF, INDEX, WORK, IWORK) -C OUTPUT U - NUMERICAL SOLUTION: +C OUTPUT U - NUMERICAL SOLUTION: DO I = 1, NP+1 - XI = XBYS(I)*DEXP(-RRT) - (1.0D0-DEXP(-RRT))/RRT + XI(1) = XBYS(I)*DEXP(-RRT) - (1.0D0-DEXP(-RRT))/RRT CALL VALUES(XI, USOL, SCRTCH, 1, 1, 1, 0, WORK) PRICEBYS(I) = USOL(1,1,1) - ENDDO - + ENDDO + RETURN END - + C ------------------------------------------------------------------------------ SUBROUTINE F(T, X, U, UX, UXX, FVAL, NPDE) - + IMPLICIT REAL*8 (A-H, O-Z) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE), FVAL(NPDE) COMMON /GEAR0/ HUSED, NQUSED, NS, NF, NJ COMMON /PARAMS/ PI COMMON /ASIAN1/ SIGMAT, RRT, XM, Z, MODEL COMMON /ASIAN2/ SIGMA, TIME, RR, XS, SS, ETA, XL, XR - + IF (MODEL.EQ.1) THEN - FR = (1.0D0-DEXP(-RR*T))/RRT + FR = (1.0D0-DEXP(-RR*T))/RRT FVAL(1) = (0.5D0*SIGMAT*SIGMAT) * ((X+FR)**2) * UXX(1) ENDIF - + IF (MODEL.EQ.2) THEN RT = (1.0D0-DEXP(-RR*T))/RR - PF = (X*SIGMA*SIGMA)/(4.0D0*DSQRT(PI*ETA)) + PF = (X*SIGMA*SIGMA)/(4.0D0*DSQRT(PI*ETA)) FVAL(1) = (0.5D0*SIGMA*SIGMA) * ((X+RT)**2) * UXX(1) FVAL(1) = FVAL(1) + PF*DEXP(-0.25D0*X*X/ETA)*(X+2.0D0*RT) ENDIF - + RETURN - + END - + C ------------------------------------------------------------------------------ - + SUBROUTINE BNDRY(T, X, U, UX, DBDU, DBDUX, DZDT, NPDE) IMPLICIT REAL*8 (A-H, O-Z) @@ -239,7 +249,7 @@ COMMON /ASIAN1/ SIGMAT, RRT, XM, Z, MODEL COMMON /ASIAN2/ SIGMA, TIME, RR, XS, SS, ETA, XL, XR -C LEFT/RIGHT BOUNDARY MODEL 1: +C LEFT/RIGHT BOUNDARY MODEL 1: IF (MODEL.EQ.1) THEN IF (X.LE.-XM) THEN DBDU (1,1) = (-Z*X + DABS(X) ) / 2.0D0 @@ -252,10 +262,10 @@ DBDUX(1,1) = 0.0D0 DZDT (1) = 0.0D0 RETURN - ENDIF + ENDIF ENDIF - -C LEFT/RIGHT BOUNDARY MODEL 2: + +C LEFT/RIGHT BOUNDARY MODEL 2: IF (MODEL.EQ.2) THEN EPS = 1.0D-20 IF (X.LE.XL ) THEN @@ -263,7 +273,7 @@ DBDUX(1,1) = 0.0D0 DZDT (1) = 0.0D0 RETURN - ENDIF + ENDIF IF (X.GE.XR ) THEN DBDU (1,1) = EPS DBDUX(1,1) = 0.0D0 @@ -271,37 +281,37 @@ RETURN ENDIF ENDIF - + RETURN - + END - + C ------------------------------------------------------------------------------ - + SUBROUTINE UINIT(X, U, NPDE) - + IMPLICIT REAL*8 (A-H, O-Z) - DIMENSION U(NPDE) + DIMENSION U(NPDE) COMMON /ASIAN1/ SIGMAT, RRT, XM, Z, MODEL COMMON /ASIAN2/ SIGMA, TIME, RR, XS, SS, ETA, XL, XR - -C NOTE : Z=+1 FOR A CALL AND Z-1 FOR A PUT - + +C NOTE : Z=+1 FOR A CALL AND Z-1 FOR A PUT + IF (MODEL.EQ.1) THEN U(1) = ( (-Z*X) + DABS(-X) ) / 2.0D0 ENDIF - + IF (MODEL.EQ.2) THEN U(1) = 0.0D0 ENDIF - + RETURN END - + C ------------------------------------------------------------------------------ - + SUBROUTINE DERIVF(T, X, U, UX, UXX, DFDU, DFDUX, DFDUXX, NPDE) - + IMPLICIT REAL*8 (A-H, O-Z) DIMENSION U(NPDE), UX(NPDE), UXX(NPDE) DIMENSION DFDU(NPDE,NPDE), DFDUX(NPDE,NPDE), DFDUXX(NPDE,NPDE) @@ -332,22 +342,22 @@ C PDE DEFINING FUNCTION F WITH RESPECT TO THE C VARIABLE UXX(J). C NOTE... THE INCOMING VALUE OF X WILL BE A COLLOCATION POINT -C VALUE. +C VALUE. PI = 4.0 * DATAN(1.0D0) - + IF (MODEL.EQ.1) THEN RT = (1.0D0-EXP(-RRT*T))/RRT DFDU(1,1) = 0.0D0 DFDUX(1,1) = 0.0D0 DFDUXX(1,1) = (SIGMAT**2) * ( X + RT ) ENDIF - + IF (MODEL.EQ.1) THEN - RT = (1.0D0-DEXP(-RR*T))/RR - F1 = (X*SIGMA*SIGMA)/(4.0D0*DSQRT(PI*ETA)) + RT = (1.0D0-DEXP(-RR*T))/RR + F1 = (X*SIGMA*SIGMA)/(4.0D0*DSQRT(PI*ETA)) F2 = DEXP(-0.25D0*X*X/ETA) - F3 = (X+2.0D0*RT) + F3 = (X+2.0D0*RT) DF1 = F1 / X DF2 = -2.0D0 * X * F2 / (4.0D0*ETA) DF3 = 1.0D0 @@ -355,23 +365,25 @@ DFDUX(1,1) = 0.0D0 DFDU(1,1) = DF1*F2*F3 + F1*DF2*F3 + F1*F2*DF3 ENDIF - + RETURN END - -C ############################################################################## - - + +C ############################################################################## + C ALGORITHM 540R (REMARK ON ALG.540), COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 18, NO. 3, SEPTEMBER, 1992, PP. 343-344. C C - SUBROUTINE PDECOL(T0, TOUT, DT, XBKPT, EPS, NINT, KORD, + SUBROUTINE PDECOL(T0, TOUT, DT, XBKPT, EPS, NINT, KORD, * NCC, NPDE, MF, INDEX, WORK, IWORK) - + IMPLICIT REAL*8 (A-H, O-Z) + + EXTERNAL REALPR, INTPR + C C C------------------------------------------------------------------------------- @@ -1171,8 +1183,10 @@ GO TO 100 C 90 DTMX = DT - 100 IF ((T+DTC) .EQ. T) WRITE(LOUT,110) - 110 FORMAT(36H WARNING.. T + DT = T ON NEXT STEP.) + 100 IF ((T+DTC) .EQ. T) CALL REALPR("WARNING.. T + DT = T ON NEXT + $ STEP.", -1, T, 0) +C 100 IF ((T+DTC) .EQ. T) WRITE(LOUT,110) +C 110 FORMAT(36H WARNING.. T + DT = T ON NEXT STEP.) C----------------------------------------------------------------------- C TAKE A TIME STEP BY CALLING THE INTEGRATOR. C----------------------------------------------------------------------- @@ -1229,63 +1243,95 @@ C TO RECOVER, DT AND DTMN ARE REDUCED BY A FACTOR OF .1 UP TO 10 C TIMES BEFORE GIVING UP. C----------------------------------------------------------------------- - 160 WRITE (LOUT,170) T - 170 FORMAT(//35H KFLAG = -1 FROM INTEGRATOR AT T = ,E16.8/ - * 41H ERROR TEST FAILED WITH DABS(DT) = DTMIN/) + 160 CALL REALPR("\n\nKFLAG = -1 FROM INTEGRATOR AT T = ", -1, T, 1) + CALL REALPR("ERROR TEST FAILED WITH DABS(DT) = DTMIN", -1, T, 0) +c$$$ 160 WRITE (LOUT,170) T +c$$$ 170 FORMAT(//35H KFLAG = -1 FROM INTEGRATOR AT T = ,E16.8/ +c$$$ * 41H ERROR TEST FAILED WITH DABS(DT) = DTMIN/) 180 IF (NHCUT .EQ. 10) GO TO 200 NHCUT = NHCUT + 1 DTMN = .1*DTMN DTC = .1*DTC - WRITE (LOUT,190) DTC - 190 FORMAT(25H DT HAS BEEN REDUCED TO ,E16.8, - * 26H AND STEP WILL BE RETRIED//) + CALL REALPR("DT HAS BEEN REDUCED TO", -1, DTC, 1) + CALL REALPR("AND STEP WILL BE RETRIED", -1, DTC, 0) +c$$$ WRITE (LOUT,190) DTC +c$$$ 190 FORMAT(25H DT HAS BEEN REDUCED TO ,E16.8, +c$$$ * 26H AND STEP WILL BE RETRIED//) JSTART = -1 GO TO 100 C - 200 WRITE (LOUT,210) - 210 FORMAT(//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//) + 200 CALL REALPR("\n\nPROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT", -1, + $ T, 0) +c$$$ 200 WRITE (LOUT,210) +c$$$ 210 FORMAT(//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//) GO TO 340 C - 220 WRITE (LOUT,230) T,DTC - 230 FORMAT(//35H KFLAG = -2 FROM INTEGRATOR AT T = ,E16.8,6H DT =, - * E16.8/52H THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//) + 220 CALL REALPR("\n\nKFLAG = -2 FROM INTEGRATOR AT T = ", -1, T, 1) + CALL REALPR("DT =", -1, DT, 1) + CALL REALPR("THE REQUESTED ERROR IS SMALLER THAN CAN BE + $ HANDLED", -1, T, 0) +c$$$ 220 WRITE (LOUT,230) T,DTC +c$$$ 230 FORMAT(//35H KFLAG = -2 FROM INTEGRATOR AT T = ,E16.8,6H DT =, +c$$$ * E16.8/52H THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//) GO TO 340 C - 240 WRITE (LOUT,250) T - 250 FORMAT(//37H INTEGRATION HALTED BY DRIVER AT T = ,E16.8/ - * 56H EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION/) + 240 CALL REALPR("\n\nINTEGRATION HALTED BY DRIVER AT T = ", -1, T, 1) + CALL REALPR("EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE + $ PRECISION", -1, T, 0) +c$$$ 240 WRITE (LOUT,250) T +c$$$ 250 FORMAT(//37H INTEGRATION HALTED BY DRIVER AT T = ,E16.8/ +c$$$ * 56H EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION/) KFLAG = -2 GO TO 340 C - 260 WRITE (LOUT,270) T - 270 FORMAT(//35H KFLAG = -3 FROM INTEGRATOR AT T = ,E16.8/ - * 45H CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED/) + 260 CALL REALPR("\n\nKFLAG = -3 FROM INTEGRATOR AT T =", -1, T, 1) + CALL REALPR("CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED", -1 , T, + $ 0) +c$$$ 260 WRITE (LOUT,270) T +c$$$ 270 FORMAT(//35H KFLAG = -3 FROM INTEGRATOR AT T = ,E16.8/ +c$$$ * 45H CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED/) GO TO 180 C - 280 WRITE (LOUT,290) - 290 FORMAT(//28H SINGULAR MATRIX ENCOUNTERED, - * 35H PROBABLY DUE TO STORAGE OVERWRITES//) + 280 CALL REALPR("\n\nSINGULAR MATRIX ENCOUNTERED", -1, T, 0) + CALL REALPR("PROBABLY DUE TO STORAGE OVERWRITES", -1, T, 0) +c$$$ 280 WRITE (LOUT,290) +c$$$ 290 FORMAT(//28H SINGULAR MATRIX ENCOUNTERED, +c$$$ * 35H PROBABLY DUE TO STORAGE OVERWRITES//) KFLAG = -4 GO TO 340 C - 300 WRITE(LOUT,310) T,TOUT,DTC - 310 FORMAT(//45H INDEX = -1 ON INPUT WITH (T-TOUT)*DT .GE. 0./ - * 4H T =,E16.8,9H TOUT =,E16.8,8H DTC =,E16.8/ - * 44H INTERPOLATION WAS DONE AS ON NORMAL RETURN./ - * 41H DESIRED PARAMETER CHANGES WERE NOT MADE.) + 300 CALL REALPR("\n\nINDEX = -1 ON INPUT WITH (T-TOUT)*DT .GE. 0." , + $ -1, T, 0) + CALL REALPR("T =", -1, T, 1) + CALL REALPR("TOUT =", -1, TOUT, 1) + CALL REALPR("DTC =", -1, DTC, 1) + CALL REALPR("INTERPOLATION WAS DONE AS ON NORMAL RETURN.", -1, T, + $ 0) + CALL REALPR("DESIRED PARAMETER CHANGES WERE NOT MADE.", -1, T, 0) +c$$$ 300 WRITE(LOUT,310) T,TOUT,DTC +c$$$ 310 FORMAT(//45H INDEX = -1 ON INPUT WITH (T-TOUT)*DT .GE. 0./ +c$$$ * 4H T =,E16.8,9H TOUT =,E16.8,8H DTC =,E16.8/ +c$$$ * 44H INTERPOLATION WAS DONE AS ON NORMAL RETURN./ +c$$$ * 41H DESIRED PARAMETER CHANGES WERE NOT MADE.) CALL INTERP(TOUT,WORK(IW10),NEQN,WORK(IW6)) INDEX = -5 RETURN C - 320 WRITE(LOUT,330) IERID - 330 FORMAT(//24H ILLEGAL INPUT...INDEX= ,I3//) + 320 CALL INTPR("\n\nILLEGAL INPUT...INDEX=", -1, IERID, 1) +c$$$ 320 WRITE(LOUT,330) IERID +c$$$ 330 FORMAT(//24H ILLEGAL INPUT...INDEX= ,I3//) INDEX = IERID RETURN C - 335 WRITE(LOUT,336) IWSTOR,IWSAVE,IISTOR,IISAVE - 336 FORMAT(//21H INSUFFICIENT STORAGE/24H WORK MUST BE OF LENGTH, - * I10,5X,12HYOU PROVIDED,I10/24H IWORK MUST BE OF LENGTH,I10,5X, - * 12HYOU PROVIDED,I10//) + 335 CALL INTPR("\n\nINSUFFICIENT STORAGE", -1, IWSTOR, 0) + CALL INTPR("WORK MUST BE OF LENGTH", -1, IWSTOR, 1) + CALL INTPR("YOU PROVIDED", -1, IWSAVE, 1) + CALL INTPR("IWORK MUST BE OF LENGTH", -1, IISTOR, 1) + CALL INTPR("YOU PROVIDED", -1, IISAVE, 1) +c$$$ 335 WRITE(LOUT,336) IWSTOR,IWSAVE,IISTOR,IISAVE +c$$$ 336 FORMAT(//21H INSUFFICIENT STORAGE/24H WORK MUST BE OF LENGTH, +c$$$ * I10,5X,12HYOU PROVIDED,I10/24H IWORK MUST BE OF LENGTH,I10,5X, +c$$$ * 12HYOU PROVIDED,I10//) INDEX = IERID RETURN C @@ -1407,7 +1453,7 @@ C*** C*** DATA LOUT,NOGAUS,MAXDER,UROUND/6, 0, 5, 5.960464D-08/ C - DATA LOUT,NOGAUS,MAXDER,UROUND/66, 0, 5, 2.22D-16/ + DATA LOUT,NOGAUS,MAXDER,UROUND/66, 0, 5, 2.22D-16/ END C C @@ -2507,7 +2553,7 @@ C C ############################################################################## C -C +C SUBROUTINE DIFFUN( N, T, Y, YDOT, IER, PW, IPIV, WORK, IWORK ) IMPLICIT REAL*8 (A-H, O-Z) C------------------------------------------------------------------------------- @@ -2545,7 +2591,7 @@ C C ############################################################################## C -C +C SUBROUTINE ADDA( PW,N0,A,ILEFT,BC,NPDE ) IMPLICIT REAL*8 (A-H, O-Z) C------------------------------------------------------------------------------- @@ -2598,7 +2644,7 @@ C C ############################################################################## C -C +C SUBROUTINE RES( T,H,C,V,R,NPDE,NCPTS,A,ILEFT,BC,DBDU,DBDUX,DZDT, * XC,UVAL ) IMPLICIT REAL*8 (A-H, O-Z) @@ -2848,7 +2894,7 @@ COMMON /OPTION/ NOGAUS,MAXDER COMMON /GEAR1/ T,H,DUMMY(4),N,IDUMMY(2),JSTART DIMENSION Y0(NEQN),Y(NEQN,MAXDER+1) - + DO 10 I = 1,N 10 Y0(I) = Y(I,1) L = JSTART + 1 @@ -3193,6 +3239,6 @@ RETURN END - + C ------------------------------------------------------------------------------