Thank you for your help, even though there was such an obvious mistake, Im sorry for that I have now tried to incorporate your suggested solution, but just as last time (the other post that you referred to), I get the values of the initial parameters when I run nlminb. I have changed the code a bit, see below, due to some error messages. Can anyone see what I am doing wrong?
Thank you in advance! http://r.789695.n4.nabble.com/file/n3758834/S%26P_500_calls%2C_jan-jun_2010.csv S%26P_500_calls%2C_jan-jun_2010.csv setwd("F:/Data til speciale/") ############## Calibration of Heston model parameters marketdata <- read.csv(file="S&P 500 calls, jan-jun 2010.csv", header=TRUE, sep=";") BS_Call <- function(S0, K, T, r, sigma, q) { sig <- sigma * sqrt(T) d1 <- (log (S0/K) + (r - q + sigma^2/2) * T ) / sig d2 <- d1 - sig Presentvalue <- exp(-r*T) return (S0 * exp(-q*T) * pnorm(d1) - K*Presentvalue*pnorm(d2)) } #------------------------- Values ---------------------------------- #### Data imported S0 <- 1136.02 X <- marketdata[1:460,9] t <- marketdata[1:460,17]/365 #Notice the T is measured in years now implvol <- marketdata[1:460,12] k <- log(X/(S0*exp(r-q)*t)) ###### Initial values kappa <- 0.0663227 # Lambda = -kappa rho <- -0.6678461 eta <- 0.002124704 theta <- 0.0001421415 v0 <- 0.0001421415 q <- 0.02145608 r <- 0.01268737 #----------------------------------------------------------------------------- #### The price of a Call option (Eq. (5.6) of The Volatility Surface, Gatheral) # In terms of log moneyness Price_call <- function(phi, k, t) { integrand <- function(u) {Re(exp(-1i*u*k)*phi(u - 1i/2, t)/(u^2 + 1/4))} res <- S0*exp(-q*t) - exp(k/2)/pi*integrate(Vectorize(integrand),lower=0,upper=Inf, subdivisions=460)$value return(res) } # The characteric formula for the Heston model (Eq. XX) phiHeston <- function(kappa, rho, eta, theta, v0) { lambda <- -kappa function(u, t) { alpha <- -u*u/2 - 1i*u/2 beta <- lambda - rho*eta*1i*u gamma <- eta^2/2 d <- sqrt(beta*beta - 4*alpha*gamma) rplus <- (beta + d)/(eta^2) rminus <- (beta - d)/(eta^2) g <- rminus / rplus D <- rminus * (1 - exp(-d*t))/ (1 - g*exp(-d*t)) C <- lambda* (rminus * t - 2/eta^2 * log( (1 - g*exp(-(d*t)))/(1 - g)) ) return(exp(C*theta + D*v0)) } } ## Calculating the Heston model price with fourier HestonCall<-function(k,t) { res<-Price_call(phiHeston(kappa, rho, eta, theta, v0),k,t) return(res) } ##### Vectorizing the function to handle vectors of strikes and maturities HestonCallVec <- function(k,t) { mapply (HestonCall, k, t) } lb <- c(0, -0.9, 0, 0, 0) ub <- c(100, 0.9, 0.5, 1, 1) #BS_Call(S0, exp(k), t, r, implvol, q) difference <- function(k, t, S0, r, implvol, q, kappa, rho, eta, theta, v0) { return(HestonCallVec(k,t) - BS_Call(S0, exp(k), t, r, implvol, q)) } difference(k,t,S0, r, implvol, q, kappa, rho, eta, theta, v0) y <- function(par) {kappahat<-par[1]; rhohat<-par[2]; etahat<-par[3]; thetahat<-par[4]; v0hat<-par[5]; sum(difference(k, t, S0, r, implvol, q, kappahat, rhohat, etahat, thetahat, v0hat)^2)} nlminb(start=c(kappa, rho, eta, theta, v0), objective = y, lower =lb, upper =ub) -- View this message in context: http://r.789695.n4.nabble.com/Error-message-object-of-type-closure-is-not-subsettable-tp3752886p3758834.html Sent from the R help mailing list archive at Nabble.com. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
