# R/ModelHeston.R In mfrdixon/MLEMVD: Maximum Likelihood Estimation for Diffusion Equations

#### Documented in ModelHeston

#' Compute the maximum likelihood estimate of the Heston model
#'
#' Heston model:
#' \eqn{dln(S_t) = (\mu - \sigma^2/2) dt + \sqrt{V_t}dW_t^{1}}
#' \eqn{dV_t = \kappa(\theta - V_t)dt  + \sigma \sqrt{V_t}dW_t^{2}}
#' 4 parameters to be estimated: (rho, kappa, theta, sigma)
#'
#' @param x Observation of the state variable at time t
#' @param x0 Observation of the state variable at time t-1
#' @param del The time step between the current and previous observation
#' @param param The parameter 4-vector (rho,kappa>0,theta>0,sigma>0)
#' @param args Optional additional arguments
#' @export output a list with a llk variable storing the result of the log likelihood calculation
#' @examples
#' \dontrun{
#' ModelHeston(0.4,0.3,0.1,c(0.1,0.3,0.2,0.1))
#' }

# The Heston price, vega and the implied volatility approximation is computed in c++.
# filename<-paste(getwd(),'/src/HestonFourierCosine.cpp',sep='')
# sourceCpp(filename)

ModelHeston<- function(x,x0,del,param,args=NULL){

output <- list()
rho   <- param[1]
kappa <- param[2]
theta <- param[3]
sigma <- param[4]

if (!is.null(args)){

if (args$mode=='implied') { S<-exp(x[1]) K<-S objfun<-function(v){ objfun<- HestonCOS(S,K,args$T_0,args$rate,args$q,param[4],param[2],param[3],v,param[1],args$callput, args$N)
}

args$v_0 <- getImpliedVolatility(S,x[3],K,args$T_0,args$rate,args$q,param[4],param[2],param[3],args$v_0,param[1],args$callput, 0.01,100, args$N) # the implied volatility #print(paste('implied vol estimate: ',v_0)) # calculate the vega to obtain the Jacobian output$v <- args$v_0 dVdv0 <- HestonVega(S,K,args$T_0,args$rate,args$q,sigma,kappa,theta,args$v_0,rho,args$callput,args$N) #print(paste('Jacobian: ',dVdv0)) #print(paste('Numerical Jacobian', grad(objfun,x[2]))) J <- dVdv0 if (is.nan(log(J))){ #print(paste('Error: NAN occured in ModelHeston.R: ',J,sigma,kappa,theta,rho,args$v_0,args$T_0)) output$llk <-0
}
else
output$llk <- -log(J) } } # m,a,b,s,r param_prime <- c(args$rate-args$q, kappa*theta, kappa, sigma, rho) model <- ModelB6(x,x0,del,param_prime)$llk
if(is.nan(model)){
print("Warning: Model likelihood is nan")
print(param_prime)
model <- 0
}
if (is.null(output$llk)) { output$llk <- 0
}

output$llk <- output$llk + model

return(output)
}
mfrdixon/MLEMVD documentation built on May 1, 2018, 11:38 p.m.