Nothing
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General
# Public License along with this library; if not, write to the
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA 02111-1307 USA
################################################################################
# FUNCTION: EBM DENSITY APPROXIMATIONS:
# dlognorm log-Normal density an derivatives
# plognorm log-Normal, synonyme for plnorm
# dgam Gamma density, synonyme for dgamma
# pgam Gamma probability, synonyme for pgamma
# drgam Reciprocal-Gamma density
# prgam Reciprocal-Gamma probability
# djohnson Johnson Type I density
# pjohnson Johnson Type I probability
# FUNCTION : MOMENTS FOR EBM DENSITY APPROXIMATIONS:
# mnorm Moments of Normal density
# mlognorm Moments of log-Normal density
# mrgam Moments of reciprocal-Gamma density
# mjohnson Moments of the Johnson Type-I density
# masian Moments of Asian Option density
# .DufresneMoments Internal Function used by masian()
# FUNCTION: NUMERICAL DERIVATIVES:
# derivative First and second numerical derivative
# FUNCTION: ASIAN DENSITY:
# d2EBM Double Integrated EBM density
# .thetaEBM Internal Function used to compute *2EBM()
# .psiEBM Internal Function used to compute *2EBM()
# dEBM Exponential Brownian motion density
# pEBM Exponential Brownian motion probability
# .gxuEBM Internal Function used to compute *EBM()
# .gxtEBM Internal Function used to compute *EBM()
# .gxtuEBM Internal Function used to compute *EBM()
# dasymEBM Exponential Brownian motion asymptotic density
################################################################################
################################################################################
# EBM DENSITY APPROXIMATIONS:
dlognorm =
function(x, meanlog = 0, sdlog = 1, deriv = c(0, 1, 2))
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the log-Normal density or its first or
# second derivative.
# Arguments:
# Details:
# Uses the function dlnorm().
# See also:
# dlnorm(x, meanlog = 0, sdlog = 1, log = FALSE)
# plnorm(q, meanlog = 0, sdlog = 1, lower.tail = TRUE, log.p = FALSE)
# qlnorm(p, meanlog = 0, sdlog = 1, lower.tail = TRUE, log.p = FALSE)
# rlnorm(n, meanlog = 0, sdlog = 1)
# FUNCTION:
# Settings:
deriv = deriv[1]
# Function:
result = dlnorm(x, meanlog = meanlog, sdlog = sdlog)
# First derivative, if desired:
if (deriv == 1) {
h1 = -(1/x + (log(x)-meanlog)/(sdlog^2*x))
result = result * h1
}
# Second derivative, if desired:
if (deriv == 2) {
h1 = -(1/x + (log(x)-meanlog)/(sdlog^2*x))
h2 = -(-1/x^2 + (-1/x^2)*(log(x)-meanlog)/sdlog^2 + 1/(sdlog^2*x^2))
result = result * (h1^2 + h2)
}
# Return Value:
result
}
# ------------------------------------------------------------------------------
plognorm =
function(q, meanlog = 0, sdlog = 1)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the log-Normal probability.
# Details:
# Uses the function plnorm().
# See also:
# dlnorm(x, meanlog = 0, sdlog = 1, log = FALSE)
# plnorm(q, meanlog = 0, sdlog = 1, lower.tail = TRUE, log.p = FALSE)
# qlnorm(p, meanlog = 0, sdlog = 1, lower.tail = TRUE, log.p = FALSE)
# rlnorm(n, meanlog = 0, sdlog = 1)
# FUNCTION:
# Resul:
result = plnorm(q = q, meanlog = meanlog, sdlog = sdlog)
# Return Value:
result
}
# ******************************************************************************
dgam =
function(x, alpha, beta)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the Gamma density.
# Details:
# The function is a synonym to "dgamma".
# See also:
# dgamma(x, shape, rate=1, scale=1/rate, log = FALSE)
# pgamma(q, shape, rate=1, scale=1/rate, lower.tail = TRUE, log = FALSE)
# qgamma(p, shape, rate=1, scale=1/rate, lower.tail = TRUE, log = FALSE)
# rgamma(n, shape, rate=1, scale=1/rate)
# FUNCTION:
# Return Value:
dgamma(x = x, shape = alpha, scale = beta)
}
# ------------------------------------------------------------------------------
pgam =
function(q, alpha, beta, lower.tail = TRUE)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the Gamma probability.
# Details:
# The function is a synonym to "pgamma".
# See also:
# dgamma(x, shape, rate=1, scale=1/rate, log = FALSE)
# pgamma(q, shape, rate=1, scale=1/rate, lower.tail = TRUE, log.p = FALSE)
# qgamma(p, shape, rate=1, scale=1/rate, lower.tail = TRUE, log.p = FALSE)
# rgamma(n, shape, rate=1, scale=1/rate)
# FUNCTION:
# Return Value:
pgamma(q = q, shape = alpha, scale = beta, lower.tail = lower.tail)
}
# ------------------------------------------------------------------------------
drgam =
function(x, alpha = 1, beta = 1, deriv = c(0, 1, 2))
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the reciprocal-Gamma density.
# See also:
# dgamma(x, shape, rate=1, scale=1/rate, log = FALSE)
# pgamma(q, shape, rate=1, scale=1/rate, lower.tail = TRUE, log.p = FALSE)
# qgamma(p, shape, rate=1, scale=1/rate, lower.tail = TRUE, log.p = FALSE)
# rgamma(n, shape, rate=1, scale=1/rate)
# FUNCTION:
# Function Value:
deriv = deriv[1]
gr = dgamma(x = 1/x, shape = alpha, scale = beta) / (x^2)
result = gr
# First Derivative:
if (deriv == 1) {
h = function(x, alpha, beta) {
-(alpha+1)/x + 1/(beta*x^2)
}
gr1 = gr*h(x, alpha, beta)
result = gr1
}
# Second Derivative:
if (deriv == 2) {
h = function(x, alpha, beta) {
-(alpha+1)/x + 1/(beta*x^2)
}
h1 = function(x, alpha, beta) {
+(alpha+1)/x^2 - 2/(beta*x^3)
}
gr2 = gr*(h(x, alpha, beta)^2 + h1(x, alpha, beta))
result = gr2
}
# Return Value:
result
}
# ------------------------------------------------------------------------------
prgam =
function(q, alpha = 1, beta = 1, lower.tail = TRUE)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the reciprocal-Gamma probability
# See also:
# dgamma(x, shape, rate=1, scale=1/rate, log = FALSE)
# pgamma(q, shape, rate=1, scale=1/rate, lower.tail = TRUE, log.p = FALSE)
# qgamma(p, shape, rate=1, scale=1/rate, lower.tail = TRUE, log.p = FALSE)
# rgamma(n, shape, rate=1, scale=1/rate)
# FUNCTION:
# Return Value:
1 - pgamma(q = 1/q, shape = alpha, scale = beta, lower.tail = lower.tail)
}
# ------------------------------------------------------------------------------
djohnson =
function(x, a = 0, b = 1, c = 0, d = 1, deriv = c(0, 1, 2))
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the Johnson Type-I density.
# FUNCTION:
# Function Value:
deriv = deriv[1]
z = a + b * log( (x-c)/d )
z1 = b / (x-c)
phi = dnorm(z, mean = 0, sd = 1)
johnson = phi * z1
result = johnson
# First Derivative:
if (deriv == 1) {
z2 = -b / (x-c)^2
johnson1 = phi * ( z2 - z*z1^2 )
result = johnson1
}
# Second Derivative:
if (deriv == 2) {
z2 = -b / (x-c)^2
z3 = 2 * b / (x-c)^3
johnson2 = phi * ( - z*z1*z2 + z^2*z1^3 + z3 -z1^3 - 2*z*z1*z2 )
result = johnson2
}
# Return Value:
result
}
# ------------------------------------------------------------------------------
pjohnson =
function(q, a = 0, b = 1, c = 0, d = 1)
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates the Johnson Type-I probability.
# FUNCTION:
# Type I:
z = a + b * log( (q-c) / d )
# Return Value:
pnorm(q = z, mean = 0, sd = 1)
}
################################################################################
# MOMENTS FOR EBM DENSITY APPROXIMATIONS:
mnorm =
function(mean = 0, sd = 1)
{ # A function implemented by Diethelm Wuertz
# Description:
# Computes the moments for the Normal distribution.
# FUNCTION:
# Raw Moments:
M = c(
mean, mean^2+sd^2, mean*(mean^2+3*sd*2), mean^4+6*mean^2*sd^2+3*sd^4 )
# Centered Moments:
m = M
m[2] = M[2] - 2*M[1]*M[1] + M[1]^2
m[3] = M[3] - 3*M[2]*M[1] + 3*M[1]*M[1]^2 - M[1]^3
m[4] = M[4] - 4*M[3]*M[1] + 6*M[2]*M[1]^2 - 4*M[1]*M[1]^3 + M[1]^4
# Fischer Parameters - Skewness and Kurtosis:
f = c(NA, NA)
f[1] = m[3] / m[2]^(3/2)
f[2] = m[4] / m[2]^2 - 3
# Return Value:
list(rawMoments = M, centralMoments = m, fisher = f)
}
# ------------------------------------------------------------------------------
mlognorm =
function(meanlog = 0, sdlog = 1)
{ # A function implemented by Diethelm Wuertz
# Description:
# Computes the moments for the Log-Normal distribution.
# FUNCTION:
# Raw Moments:
n = 1:4
M = exp ( n * meanlog + n^2 * sdlog^2/2 )
# Centered Moments:
m = M
m[2] = M[2] - 2*M[1]*M[1] + M[1]^2
m[3] = M[3] - 3*M[2]*M[1] + 3*M[1]*M[1]^2 - M[1]^3
m[4] = M[4] - 4*M[3]*M[1] + 6*M[2]*M[1]^2 - 4*M[1]*M[1]^3 + M[1]^4
# Fischer Parameters - Skewness and Kurtosis:
f = c(NA, NA)
f[1] = m[3] / m[2]^(3/2)
f[2] = m[4] / m[2]^2 - 3
# Return Value:
list(rawMoments = M, centralMoments = m, fisher = f)
}
# ------------------------------------------------------------------------------
mrgam =
function(alpha = 1/2, beta = 1)
{ # A function implemented by Diethelm Wuertz
# Description:
# Computes the moments for the Reciprocal-Gamma distribution.
# FUNCTION:
# Raw Moments:
M = rep(0, times = 4)
M[1] = 1 / (beta*(alpha - 1))
M[2] = M[1] / (beta*(alpha - 2))
M[3] = M[2] / (beta*(alpha - 3))
M[4] = M[3] / (beta*(alpha - 4))
# Centered Moments:
m = M
m[2] = M[2] - 2*M[1]*M[1] + M[1]^2
m[3] = M[3] - 3*M[2]*M[1] + 3*M[1]*M[1]^2 - M[1]^3
m[4] = M[4] - 4*M[3]*M[1] + 6*M[2]*M[1]^2 - 4*M[1]*M[1]^3 + M[1]^4
# Fischer Parameters - Skewness and Kurtosis:
f = c(NA, NA)
f[1] = m[3] / m[2]^(3/2)
f[2] = m[4] / m[2]^2 - 3
# Return Value:
list(rawMoments = M, centralMoments = m, fisher = f)
}
# ------------------------------------------------------------------------------
mjohnson =
function(a, b, c, d)
{ # A function implemented by Diethelm Wuertz
# Description:
# Computes the moments for the Johnson Type-I distribution
# FUNCTION:
# Raw Moments:
M = c(NA, NA, NA, NA)
# Centered Moments:
m = M
m[2] = M[2] - 2*M[1]*M[1] + M[1]^2
m[3] = M[3] - 3*M[2]*M[1] + 3*M[1]*M[1]^2 - M[1]^3
m[4] = M[4] - 4*M[3]*M[1] + 6*M[2]*M[1]^2 - 4*M[1]*M[1]^3 + M[1]^4
# Fischer Parameters - Skewness and Kurtosis:
f = c(NA, NA)
f[1] = m[3] / m[2]^(3/2)
f[2] = m[4] / m[2]^2 - 3
# Return Value:
list(rawMoments = M, centralMoments = m, fisher = f)
}
# ------------------------------------------------------------------------------
masian =
function(Time = 1, r = 0.045, sigma = 0.30)
{ # A function implemented by Diethelm Wuertz
# Description:
# Computes the moments for the Asian-Option distribution
# FUNCTION:
# Raw Moments:
M = .DufresneMoments(M = 4, Time = Time, r = r, sigma = sigma)
# Centered Moments:
m = M
m[2] = M[2] - 2*M[1]*M[1] + M[1]^2
m[3] = M[3] - 3*M[2]*M[1] + 3*M[1]*M[1]^2 - M[1]^3
m[4] = M[4] - 4*M[3]*M[1] + 6*M[2]*M[1]^2 - 4*M[1]*M[1]^3 + M[1]^4
# Fischer Parameters - Skewness and Kurtosis:
f = c(NA, NA)
f[1] = m[3] / m[2]^(3/2)
f[2] = m[4] / m[2]^2 - 3
# Return Value:
list(rawMoments = M, centralMoments = m, fisher = f)
}
# ------------------------------------------------------------------------------
.DufresneMoments =
function (M = 4, Time = 1, r = 0.045, sigma = 0.30)
{ # A function implemented by Diethelm Wuertz
# Description:
# Computes raw moments for the Asian-Option distribution.
# Note:
# Called by function masian()
# FUNCTION:
# Internal Function:
moments =
function (M, tau, nu) {
d = function(j, n, beta) {
d = 2^n
for (i in 0:n) if (i != j) d = d / ( (beta+j)^2 - (beta+i)^2 )
d
}
moments = rep(0, length = M)
for (n in 1:M) {
moments[n] = 0
for (j in 0:n) moments[n] = moments[n] +
d(j, n, nu/2)*exp(2*(j^2+j*nu)*tau)
moments[n] = prod(1:n) * moments[n] / (2^(2*n))
}
moments
}
# Compute:
tau = sigma^2*Time/4
nu = 2*r/sigma^2-1
ans = (4/sigma^2)^(1:M) * moments(M, tau, nu)
# Return Value:
ans
}
################################################################################
derivative =
function(x, y, deriv = c(1, 2))
{ # A function implemented by Diethelm Wuertz
# Description:
# Calculates numerically the first or second derivative
# of the functuin y(x) by finite differences.
# Arguments:
# x - a numeric vector of values
# y - a numeric vectror of function values y(x)
# deriv - the degree of differentiation, either 1 or 2.
# FUNCTION:
# Stop in the case of wrong argument deriv:
dseriv = deriv[1]
if (deriv < 1 || deriv > 2) stop("argument error")
# Function to calculate the next derivative by differences:
"calcderiv" = function(x,y) {
list(x = x[2:length(x)]-diff(x)/2, y = diff(y)/diff(x))}
# First Numerical Derivative:
result = calcderiv(x, y)
# Second Numerical Derivative, if desired:
if (deriv == 2) result = calcderiv(result$x, result$y)
# Return Value:
list(x = result$x, y = result$y)
}
################################################################################
d2EBM =
function(u, t = 1)
{ # A function written by Diethelm Wuertz
# Description:
# Calculate the density integral "f_A_t(u)" given by
# equation 4.36 in: R. Gould, "The Distribution of the
# Integral of Exponential Brownian Motion".
# Arguments:
# t - numeric value
# u - numeric value
# FUNCTION:
# Function to be integrated:
f = function(x, tt, uu) {
fx = rep(0, length=length(x))
for (i in 1:length(x) )
fx[i] = (1/uu) * exp(-(1+exp(2*x[i]))/(2*uu)) *
.thetaEBM(r=exp(x[i])/uu, u=tt)
fx }
# Integrate:
result = rep(0, length = length(u))
for (i in 1:length(u)) {
result[i] = integrate(f, lower = -16, upper = 4, tt = t, uu = u[i],
subdivisions = 100, rel.tol=.Machine$double.eps^0.25,
abs.tol=.Machine$double.eps^0.25)$value }
# Return Value:
result
}
# ------------------------------------------------------------------------------
.thetaEBM =
function(r, u)
{ # A function written by Diethelm Wuertz
# Description:
# Calculate the integral "\theta_r(u)" given by equations
# 2.22 and 2.23 in: R. Gould, "The Distribution of the
# Integral of Exponential Brownian Motion".
# Arguments:
# r - vector of numeric values
# u - numeric value
# FUNCTION:
# Function to be integrated:
f = function(x, rr, uu) {
fx = rep(0, length = length(x))
for (i in 1:length(x) )
fx[i] = exp(-x[i]^2/(2*uu)) * exp(-rr*cosh(x[i])) *
sinh(x[i]) * sin(pi*x[i]/uu)
fx }
# Loop over r-Vector:
result = rep(0, length=length(r))
for ( i in 1: length(r) ) {
result[i] = integrate(f, lower = 0, upper = 30, rr = r[i], uu = u,
subdivisions = 100, rel.tol = .Machine$double.eps^0.25,
abs.tol = .Machine$double.eps^0.25)$value
result[i] = result[i] * (r/sqrt((2*u*pi^3))) * exp(pi^2/(2*u)) }
# Return Value:
result
}
# ------------------------------------------------------------------------------
.psiEBM =
function(r, u)
{ # A function written by Diethelm Wuertz
# Description:
# Calculate the integral "\psi_r(u)" given by equations
# 2.22 and 2.23 in: R. Gould, "The Distribution of the
# Integral of Exponential Brownian Motion".
# Arguments:
# r - vector of numeric values
# u - numeric value
# FUNCTION:
# Calculate psi() from theta():
result = sqrt(2*u*pi^3) * exp(-pi^2/(2*u)) * .thetaEBM(r, u)
# Return Value:
result
}
# ------------------------------------------------------------------------------
dEBM =
function(u, t = 1)
{ # A function written by Diethelm Wuertz
# Arguments;
# t - a numeric value
# u - a vector of numeric values
# FUNCTION:
# Calculate Density:
result = rep(0, times = length(u))
for (i in 1:length(u) ) {
result[i] = integrate(.gxtuEBM, lower = 0, upper = 100,
t = t, u = u[i])$value
}
# Return Value:
result
}
# ------------------------------------------------------------------------------
pEBM =
function(u, t = 1)
{ # A function written by Diethelm Wuertz
# Arguments;
# t - a numeric value
# u - a vector of numeric value
# FUNCTION:
# Calculate Probability:
result = rep(0, times = length(u))
result[1] = integrate(dEBM, lower = 0, upper = u[1], t = t)$value
if (length(u) > 1) {
for (i in 2:length(u) ) {
result[i] = result[i-1] + integrate(
dEBM, lower = u[i-1], upper = u[i], t = t)$value
}
}
# Return Value:
result
}
# ------------------------------------------------------------------------------
.gxuEBM =
function(x, u)
{ # A function written by Diethelm Wuertz
# Description:
# Interchange the Integrals - and first integrate:
# 1/u^2 * sinh(x) * exp(-(1/(2*u))*(1+exp(2*x))) * exp(x) *
# exp(-exp(x)*cosh(y)/u)
# FUNCTION:
# Compute g(x, u):
fx = rep(0, length = length(x))
if (u > 0) {
for ( i in 1:length(x) ) {
su = (u)^(-3/2) * sinh(x[i])
cx = cosh(x[i])/sqrt(u)
sx = sinh(x[i])/sqrt(u)
Asymptotics = exp(x[i]) / sqrt(u) / 2
if (Asymptotics <= 33.0) {
fx[i] = su * pnorm(-cx) / dnorm(sx) }
else {
fx[i] = su * exp(-1/(2*u)) * (1-1/cx^2+3/cx^4) / cx } } }
# Return Value:
fx
}
# ------------------------------------------------------------------------------
.gxtEBM =
function(x, t)
{ # A function written by Diethelm Wuertz
# Description:
# FUNCTION:
# Compute g(x, t):
fx = exp(pi^2/(2*t)) / sqrt(2*t*pi^3) * exp(-x^2/(2*t))
# Return Value:
fx
}
# ------------------------------------------------------------------------------
.gxtuEBM =
function(x, t, u)
{ # A function written by Diethelm Wuertz
# FUNCTION:
# Result:
fx = .gxtEBM(x = x, t = t) * .gxuEBM(x = x, u = u) * sin(pi*x/t)
# Return Value:
fx
}
# ------------------------------------------------------------------------------
dasymEBM =
function(u, t = 1)
{ # A function written by Diethelm Wuertz
# Description:
# Calculates the asymptotic behavior of the density
# function f of the exponential Brownian maotion
# FUNCTION:
# Asymptotic Density:
alpha = log ( 8*u*exp(-2*t) )
beta = exp ( -((log(alpha/(4*t)))^2)/(8*t) )
f = sqrt(t) * exp(t/2) * exp(-alpha^2/(8*t)) * beta
# Take care of gamma function ...
warn = options()$warn
options(warn = -1)
# f = f / (u * sqrt(u) * alpha * gamma(alpha/(4 * t)))
f = f / (u * sqrt(u) * (4*t) * gamma(alpha/(4 * t)+1))
f[is.na(f)] = 0
options(warn = warn)
# Return Value:
f
}
################################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.