phi_heston <- function(a, v0, v_t, d){
gamma_a <- sqrt(k^2 - 2 * sigma^2 * 1i*a)
gammadt <- gamma_a * (tau-t)
sqrtv0vt <- sqrt(v0*v_t)
part1 <- (gamma_a * exp(-(gamma_a - k)/2 * (tau-t)) *
(1 - exp(-k * (tau-t))))/ (k * (1- exp(- gammadt)))
part2 <- exp((v0+v_t)/(sigma^2) *
((k * (1 + exp(-k*(tau-t))))/(1-exp(-k*(tau-t))) -
(gamma_a * (1 + exp(- gammadt)))/(1-exp(- gammadt))))
part3 <- Bessel::BesselI(z = ((4 * gamma_a * sqrtv0vt)/(sigma^2) *
exp(- gammadt/2)/
(1 - exp(- gammadt))), nu = 0.5*d - 1) /
Bessel::BesselI(z = ((4 * k * sqrtv0vt)/(sigma^2) * (exp(-k*(tau-t)/2))/
(1-exp(-k*(tau-t)))), nu = 0.5*d - 1)
result <- part1 * part2 * part3
return (result)
}
integrate_gauss_laguerre <- function(f, M=5){
points_weights <- gaussquad::laguerre.quadrature.rules(M)[[M]]
int <- 0
for(i in 1:M){
int <- int + points_weights[i, 'w'] *
f(points_weights[i, 'x']) * exp(points_weights[i, 'x'])
}
return (int)
}
intv <- function(n, cf, v_t){
integrand <- function(x, phi = cf){
f2 <- function(u){
# Re(phi(u)) * sin(u * x) /u
Im(phi(u) * exp(-1i * u * x)) /u
}
return(f2)
}
## integrate to "cdf"
F_x <- function (x) {
y <- 0.5 - 1/pi * integrate_gauss_laguerre(integrand(x))
return (y)
}
## endsign
endsign <- function(f, sign = 1) {
b <- sign
while (sign * f(b) < 0) b <- 10 * b
return(b)
}
## inversion
spdf.upper =Inf
spdf.lower = -Inf
invcdf <- function(u) {
subcdf <- function(t) F_x(t) - u
if (spdf.lower == -Inf)
spdf.lower <- endsign(subcdf, -1)
if (spdf.upper == Inf)
spdf.upper <- endsign(subcdf)
# browser()
return(uniroot(subcdf, lower=spdf.lower, upper=spdf.upper)$root)
}
U <- stats::runif(n)
sapply(U, invcdf)
}
hestonea_mod <- function(S, X, r, v, theta, rho, k, sigma, t = 0, tau = 1){
ST <- NULL
d1 <- (4 * k * theta)/(sigma)^2
c0 <- (sigma^2 * (1 - exp(-k*(tau-t))))/(4*k)
# sampling V
lambda <- (4*k*exp(-k*(tau-t))*v)/(sigma^2 * (1-exp(-k*(tau-t))))
vt <- c0 * stats::rchisq(n = 1, df = d1, ncp = lambda)
# Sampling int{V}
phi <- function(a, v0=v, v_t=vt, d=d1){phi_heston(a, v0=v, v_t=vt, d=d1)}
int_v <- intv(1, cf = phi, v_t=vt)
# Sampling int{v}dw
int_vdw <- (1/sigma) * (vt - v - k * theta * (tau-t) + k * int_v)
# Sampling S
if( int_v >= 0){
m <- log(S) + (r * (tau - t) - (1/2) * int_v + rho * int_vdw)
std <- sqrt((1 - rho^2)) * sqrt(int_v)
S <- exp(m + std * rnorm(1))
v <- vt
ST <- rbind(ST,S)
} else {
v <- vt
ST <- rbind(ST,S)}
ST <- as.matrix(ST, ncol=1)
Result <- ST[nrow(ST),] - X
Result[Result <= 0] = 0
call = exp(-r*tau)*Result
lista = list('call' = call, 'Result' = Result, 'Spot' = ST)
return(lista)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.