[R] nlminb - how to avoid evaluating initial parameters infinite in integrate

[Newbie]

Dear R-users.

I am faced with a problem I dont know how to solve. 
I need to calibrate the Heston stochastic volatility model, and have (to my
own belief) created a code for calculating the prices of options by this
model. However, when I calibrate the model using NLMINB I also evaluate my
initial parameters to infinity by the integrate function, and this is wrong!
I believe that this is the reason why nlminb keeps spitting out the initial
parameters as the estimates. For 0 the values of C and D, see below, should
be 0, so that the phi is equal to 1 (which makes the parameters NOT being
evaluated by infinity). It is a bit difficult to explain, and I hope you can
understand what I want to do. 
Is it possible to incorporate a statement in the integrate() saying that at
0 the conditions are different? 
Or should this be done in some sort of loop? Please help me. 

Also, I am very sorry that this message is a "repost" on my part. In my
original message the problem was a wrong use of nlminb. I hope that this new
header will make it easier for people who have knowledge in this field to
find my post. I have searched, but have not found a way to solve this
problem - I hope you are able to help me. (the problem is in the integrate
in the Price_call function where integrand refers to the phi function)


THANK you
Rikke

my code is: 

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*exp((r-q)*T)/K) + (sigma^2/2) * T ) / sig
	d2 <- d1 - sig
	Presentvalue <- exp(-r*T)
	return (Presentvalue*(S0 * exp((r-q)*T) * pnorm(d1) - K*pnorm(d2)))
}


#------------------------- Values ----------------------------------
#### Data imported
S0 <- 1136.03
X <- marketdata[1:460,9]
t <- marketdata[1:460,17]/365		#Notice the T is measured in years now
implvol <- marketdata[1:460,12]

###### 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

smallk <- log(X/(S0*exp(r-q)*t))
parameters0 <- c(kappa, rho, eta, theta, v0)
#-----------------------------------------------------------------------------


#### The price of a Call option (Eq. (5.6) of The Volatility Surface,
Gatheral)
# In terms of log moneyness

Price_call <- function(phi, smallk, t)
{
	integrand <-  function(u) {Re(exp(-1i*u*smallk)*phi(u - 1i/2, t)/(u^2 +
1/4))}
	res <- S0*exp(-q*t) -
exp(smallk/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(parameters)
{	
	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(smallk, t)
{
res<-Price_call(phiHeston(parameters),smallk,t)
return(res)
}

##### Vectorizing the function to handle vectors of strikes and maturities
HestonCallVec <- function(smallk,t)
{
mapply (HestonCall, smallk, t)
}


lb <- c(0, -1, 0, 0, 0)
ub <- c(Inf, 1, Inf, Inf, 2)


difference <- function(smallk, t, S0, r, implvol, q, parameters)
{
return(HestonCallVec(smallk,t) - BS_Call(S0, exp(smallk), t, r, implvol, q))
}

y <- function(x) {kappa<-x[1]; rho<-x[2]; eta<- x[3]; theta<- x[4];
v0<-x[5]; sum(difference(smallk, t, S0, r, implvol, q, x)^2/BS_Call(S0,
exp(smallk), t, r, implvol, q))}
nlminb(start=parameters0, objective = y, lower =lb, upper =ub)
http://r.789695.n4.nabble.com/file/n3765932/S%26P_500_calls%2C_jan-jun_2010.csv
S%26P_500_calls%2C_jan-jun_2010.csv 

--
View this message in context: http://r.789695.n4.nabble.com/nlminb-how-to-avoid-evaluating-initial-parameters-infinite-in-integrate-tp3765932p3765932.html
Sent from the R help mailing list archive at Nabble.com.