[R] How to vectorize this function

[Lynette Chang]

Here are the code and data file. I’m not sure if I put too much unrelated information here. 

My goal is to factor out volatilities from the data. I hope I can get sigV <- impVolC(callM, K, T, F, r), which has five vectors as input, and one vector as output. The length of all those six vectors are the same. However, I got stuck in the nested if-else sentence, as if-condition cannot handle vectors. I rewrite it as you guys suggestions, however, I still have one layer of if-condition. Any thoughts to improve it? Thanks.

Lynette

-------------- next part --------------


df <- read.csv(file = "/S&P500_ETF_Option_0917.csv", header = TRUE, 
               colClasses = c("integer", "character", "numeric", "numeric", "numeric",
                              "character", "numeric", "numeric", "numeric"))

call <- (df$callBid + df$callAsk)/2
put <- (df$putBid + df$putAsk)/2

y <- call - put
A <- cbind(rep(1, dim(df)[1]), -df$Strike)

x <- solve(t(A)%*%A)%*%t(A)%*%y
PVF <- x[1]
disc <- x[2]

S <- 2381



library(timeDate)
# Lets work in your environment:
getRmetricsOptions("myFinCenter")
setRmetricsOptions(myFinCenter = "America/New_York")

# define a sequence of days with  timeSequence 
t1 <- timeSequence(from = "2017-03-16", to = "2017-09-29")

# Define a calendar for your exchange (use an available one as a template, e.g. holidayNYSE) 
# subindex the sequence with isBizday using your calendar as an argument 
holidayNYSE(2017)
isBizday(t1, holidayNYSE())
t2 <- t1[isBizday(t1, holidayNYSE(2017))]

T <- length(t2)/252
q_m <- -log(PVF/S)/T
r_m <- log(disc)/(-T)


polya <- function(x){
  1/2 + sign(x)/2* sqrt(1- exp(-2*x^2/pi))
}
impVolC <- function(callM, K, T, F, r){
  
  y <- log(F/K)
  alpha <- callM/(K*exp(-r*T))
  R <- 2*alpha - exp(y) + 1
  
  A <- (exp((1 - 2/pi)*y) - exp(-(1 - 2/pi)*y))^2
  B <- 4*(exp(2/pi*y) + exp(-2/pi*y)) - 2*exp(-y)*(exp((1-2/pi)*y)+exp(-(1-2/pi)*y))*(exp(2*y) + 1 - R^2)
  C <- exp(-2*y)*(R^2 - (exp(y) -1)^2)*((exp(y) + 1)^2 - R^2)                                                
  
  beta <- (2*C)/(B + sqrt(B^2 + 4*A*C))
  gamma <- -pi/2*log(beta)
  
  if(y >= 0){
    call0 <- K*exp(-r*T)*(exp(y)*polya(sqrt(2*y)) - 0.5)
    sig <- (sqrt(gamma + y)  + ifelse(callM <= call0, -1, 1) * sqrt(gamma - y))/sqrt(T)
  }else{
    call0 <- K*exp(-r*T)*(exp(y)/2 - polya(-sqrt(-2*y)))
    sig <-  (ifelse(callM <= call0, -1, 1)*sqrt(gamma + y)  +  sqrt(gamma - y))/sqrt(T)
  }
  sig
}

F <- PVF*exp(r_m*T)
sigV <- rep(0, length(call))

for(i in 1:length(call)){
  sigV[i] <- impVolC(callM = call[i], K = df$Strike[i], T = T, F = F, r = r_m)  
}


> On Sep 20, 2018, at 1:56 PM, MacQueen, Don <macqueen1 using llnl.gov> wrote:
> 
> In addition to what the other said, if callM is a vector then an expression of the form
>   if (callM <= call0)
> is inappropriate. Objects inside the parentheses of   if()  should have length one. For example,
> 
>> if (1:5 < 3) 'a' else 'b'
> [1] "a"
> Warning message:
> In if (1:5 < 3) "a" else "b" :
>  the condition has length > 1 and only the first element will be used
> 
> 
> instead of what you have:
>         if(callM <= call0){
>           sig <- 1/sqrt(T)*(sqrt(gamma + y) - sqrt(gamma - y))
>         }else{
>           sig <- 1/sqrt(T)*(sqrt(gamma + y) + sqrt(gamma - y))
>        }
> 
> Here are a couple of (untested) possibilities:
> 
>  M.gt.0 <- callM > call0
>  sig <- 1/sqrt(T)*(sqrt(gamma + y) - sqrt(gamma - y))
>  sig[M.gt.0] <- (1/sqrt(T)*(sqrt(gamma + y) + sqrt(gamma - y)))[M.gt.0]
> 
> or
> 
>    sig <- 1/sqrt(T)*(sqrt(gamma + y)  + ifelse(callM <= call0, -1, 1) * sqrt(gamma - y))
> 
> incidentally, I would write
>   sig <- (sqrt(gamma + y) - sqrt(gamma - y))/sqrt(T)
> instead of
>   sig <- 1/sqrt(T)*(sqrt(gamma + y) - sqrt(gamma - y))
> 
> --
> Don MacQueen
> Lawrence Livermore National Laboratory
> 7000 East Ave., L-627
> Livermore, CA 94550
> 925-423-1062
> Lab cell 925-724-7509
> 
> 
> 
> On 9/20/18, 8:08 AM, "R-help on behalf of Lynette Chang" <r-help-bounces using r-project.org on behalf of momtoomax using gmail.com> wrote:
> 
>    Hello everyone,
> 
>         I’ve a function with five input argument and one output number. 
>    	  impVolC <- function(callM, K, T, F, r)
> 
>         I hope this function can take five vectors as input, then return one vector as output. My vectorization ran into problems with the nested if-else operation. As a result, I have to write another for loop to call this function. Can anyone suggest some methods to overcome it? I put my code below, thanks.
> 
>    impVolC <- function(callM, K, T, F, r){
> 
> 
>     if(y >= 0){
>         call0 <- K*exp(-r*T)*(exp(y)*polya(sqrt(2*y)) - 0.5)
>         if(callM <= call0){
>           sig <- 1/sqrt(T)*(sqrt(gamma + y) - sqrt(gamma - y))
>         }else{
>           sig <- 1/sqrt(T)*(sqrt(gamma + y) + sqrt(gamma - y))
>         }
>     }else{
>         call0 <- K*exp(-r*T)*(exp(y)/2 - polya(-sqrt(-2*y)))
>         if(callM <= call0){
>           sig <- 1/sqrt(T)*(-sqrt(gamma + y) + sqrt(gamma - y))
>         }else{
>           sig <- 1/sqrt(T)*(sqrt(gamma + y) + sqrt(gamma - y))
>         }
>     }
>     sig
>    } 
> 
>    for(i in 1:length(call)){
>     sigV[i] <- impVolC(callM = call[i], K = df$Strike[i], T = T, F = F, r = r_m)  
>    }
> 
>    ______________________________________________
>    R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>    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.
> 
>