[R] How to vectorize a function to handle two vectors

[David Winsemius]

On Aug 26, 2011, at 2:43 PM, Newbie wrote:

> Dear R-users
>
> I am trying to "vectorize" a function so that it can handle two  
> vectors of
> inputs. I want the function to use phi (a function), k[1] and t[1]  
> for the
> first price, and so on for the second, third and fourth price.

? mapply


> I tried to do
> the mapply,

Oh, you thought of that.

> but I dont know how to specify to R what input I want to be
> vectors (k and t)(see in the bottom what I tried). I have read the  
> help
> file, but I still dont understand how to do it properly. Also, I've  
> tried to
> use sapply (which seems totally wrong). However, the function uses  
> k[1] for
> t 1 to 4, and thereby returns 16 different values instead of just 4.

Huh?

> Can anyone tell me how to do this - I know the answer is simple, but  
> I dont
> understand how

I think the problem is (at least) slightly more complex than you first  
described. You should spend more time describing in unambiguous  
language the conditions and operations. Keep reading.

>
> Thank you for your time
> Kinds Rikke
>
> #------ Characteristic function of the Heston model -----#
> phiHeston <- function(parameters)
> {
> 	lambda <- parameters[1];
>  	rho <- parameters[2];
>   	eta <- parameters[3];
>  	theta <- parameters[4];
> 	v0 <- parameters[5];
>
> 	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)/(2*gamma);
>    	rminus <- (beta - d)/(2*gamma);
>    	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));
> 	}
> }
>
>
> Price_call <- function(phi, k, t)
> {
> 	integrand <-  function(u){Re(exp(-1i*u*k)*phi(u - 1i/2, t)/(u^2 +  
> 1/4))};
>    	res <- 1 - exp(k/2)/pi*integrate(integrand,lower=0,upper=Inf) 
> $value;
>    	return(res);
> }
>
> subHeston <- c(0.6067,-0.7571,0.2928,0.0707,0.0654);
> kV <- c(0.9,1,1.2,1.3)
> tV <- c(0.1,0.4,0.5, 1)
>
> HestonCallVec <- function(phi, kVec, t)
> {
> sapply(kVec, function(k){Price_call(phi, k, t)})
> }
> HestonCallVec(phiHeston(subHeston), kV, 1)
>
> subHeston <- c(0.6067,-0.7571,0.2928,0.0707,0.0654);
> kV <- c(0.9,1,1.2,1.3)
> tV <- c(0.1,0.4,0.5, 1)
>
> HestonCallVec2 <- function(phi, kVec, tVec)
> {
> sapply(tVec, function(t){HestonCallVec(phi, kVec, t)})
> }
> HestonCallVec2(phiHeston(subHeston), kV, tV)
> # This should give 4 values

Hmmmm, I get 16 values. Maybe you should put some more time into  
understanding why it does not return the structure you expect.  
phiHubHeston returns a function. Was that what you designed?

>
> This is what I tried for the mapply function, which returns a list,  
> instead
> of values.
>
> HestonCallVec <- function(phi, kVec, t)
> {
> mapply(function(k, t){Price_call(phi, k, t)})
> }
> HestonCallVec(phiHeston(subHeston), kV, tV)
> # This should give 4 values

This is what I get with the mapply call I would have expected for use  
of a function which is what phiHeston(subHeston) returns:

 > # This should give 4 values
 > mapply(phiHeston(subHeston), kV, tV)
[1] 0.9973156-0.0029158i 0.9862950-0.0124475i 0.9752339-0.0178313i
[4] 0.9406358-0.0336025i

Whether these are the "right" 4 values is for you to determine.

-- 

David Winsemius, MD
West Hartford, CT