How to loop optim over a dataframe in R

Question

I am struggling in looping the optimisation function optim to my whole dataframe. Below you can find a pice of my code where my optimisation problem works when I use scalar as Strike09 = 1142.757 but now I would like to solve this optimisation problem over my whole dataset where Strike09, Strike095, ..., blackscholesc11 are the columns of the dataset and I would like to get the 5 parameter optimising my function for each row of the dataset (corresponding to a date).

Below the function to be optimised and an example of the working optim problem when using scalar as hard parameter

error_vector_price <- function(Strike09, Strike095, Strike1, Strike105, Strike11, Liborrate, par_x, TimetoMaturity, SPX, blackscholesc09, blackscholesc095, blackscholesc1, blackscholesc105, blackscholesc11){
  
MixtureCall <- function(Strike, Liborrate, par_x, TimetoMaturity, SPX){
  
  alpha1 = log(SPX) + (par_x[2]-0.5*(par_x[4]^2)*(par_x[4]^2))*TimetoMaturity
  alpha2 = log(SPX) + (par_x[3]-0.5*(par_x[5]^2)*(par_x[5]^2))*TimetoMaturity
  beta1 = (par_x[4]^2)*(TimetoMaturity^0.5)
  beta2 = (par_x[5]^2)*(TimetoMaturity^0.5)
  theta = (par_x[1]^2)/(1+par_x[1]^2)
  
  D1 = (-log(Strike) + alpha1 + (beta1 ^ 2)) / beta1
  D2 = D1 - beta1
  D3 = (-log(Strike) + alpha2 + (beta2 ^ 2)) / beta2
  D4 = D3 - beta2
  nD1 = sapply(D1,  function(x) pnorm(x))
  nD2 = sapply(D2,  function(x) pnorm(x))
  nD3 = sapply(D3,  function(x) pnorm(x))
  nD4 = sapply(D4,  function(x) pnorm(x))
  
  call1 = exp(alpha1 + beta1 * beta1 / 2) * nD1 - Strike * nD2
  call2 = exp(alpha2 + beta2 * beta2 / 2) * nD3 - Strike * nD4
    
  InstRate = (log(1 + Liborrate * TimetoMaturity)) / TimetoMaturity
    
  Mixturecall = exp(-InstRate * TimetoMaturity) * (theta * call1 + (1 - theta) * call2)
  
  return(Mixturecall)
}

MixturePut <- function(Strike, Liborrate, par_x, TimetoMaturity, SPX){
  
  alpha1 = log(SPX) + (par_x[2]-0.5*(par_x[4]^2)*(par_x[4]^2))*TimetoMaturity
  alpha2 = log(SPX) + (par_x[3]-0.5*(par_x[5]^2)*(par_x[5]^2))*TimetoMaturity
  beta1 = (par_x[4]^2)*(TimetoMaturity^0.5)
  beta2 = (par_x[5]^2)*(TimetoMaturity^0.5)
  theta = (par_x[1]^2)/(1+par_x[1]^2)
  
  D1 = (-log(Strike) + alpha1 + (beta1 ^ 2)) / beta1
  D2 = D1 - beta1
  D3 = (-log(Strike) + alpha2 + (beta2 ^ 2)) / beta2
  D4 = D3 - beta2
  nD1 = sapply(D1, function(x) pnorm(x))
  nD2 = sapply(D2, function(x) pnorm(x))
  nD3 = sapply(D3, function(x) pnorm(x))
  nD4 = sapply(D4, function(x) pnorm(x))
  
  call1 = exp(alpha1 + beta1 * beta1 / 2) * nD1 - Strike * nD2
  call2 = exp(alpha2 + beta2 * beta2 / 2) * nD3 - Strike * nD4
    
  InstRate = (log(1 + Liborrate * TimetoMaturity)) / TimetoMaturity
    
  Mixturecall = exp(-InstRate * TimetoMaturity) * (theta * call1 + (1 - theta) * call2)
  
  fwdprice = theta * exp(alpha1 + beta1 * beta1 / 2) + (1 - theta) * exp(alpha2 + beta2 * beta2 / 2)
  
  Mixtureput = Mixturecall + (Strike - fwdprice) / (1 + Liborrate * TimetoMaturity)
  
  return(Mixtureput)
}

model_price_vector09 <- MixturePut(Strike09, Liborrate, par_x, TimetoMaturity, SPX)
model_price_vector095 <- MixturePut(Strike095, Liborrate, par_x, TimetoMaturity, SPX)
model_price_vector1 <- MixturePut(Strike1, Liborrate, par_x, TimetoMaturity, SPX)
model_price_vector105 <- MixtureCall(Strike105, Liborrate, par_x, TimetoMaturity, SPX)
model_price_vector11 <- MixtureCall(Strike11, Liborrate, par_x, TimetoMaturity, SPX)

error_vector_price <- (((blackscholesc09 - model_price_vector09)^2+(blackscholesc095 - model_price_vector095)^2+(blackscholesc1 - model_price_vector1)^2+(blackscholesc105 - model_price_vector105)^2+(blackscholesc11 - model_price_vector11)^2))

return(error_vector_price)
}

x_seed <- c(0.64,-0.57,0.25,0.53,0.36)

results <- optim( x_seed , fn = error_vector_price, gr = NULL,  Strike09 = 1142.757, Strike095 = 1206.2435, Strike1 = 1269.73, Strike105 = 1333.216, Strike11 = 1396.703, Liborrate = 0.0505750, TimetoMaturity = 0.25, SPX = 1269.73, blackscholesc09=24.995126, blackscholesc095=37.78425, blackscholesc1=57.87691, blackscholesc105=33.47135, blackscholesc11=12.671979, method = "L-BFGS-B", upper = (1))

x_star <- results$par
x_star

I now would like to loop this optimisation problem on my dataset and getting for each row the 5 parameter minimising my function. My dataset is composed by 152 obs of 13 variables, the 13 variable are the hard parameters of my optimisation problem (Strike09, Strike095, Strike1, Strike105, Strike11, Liborrate, TimetoMaturity, SPX, blackscholesc09, blackscholesc095, blackscholesc1, blackscholesc105, blackscholesc11). The idea is showed below but the loop that follows is not working and I am not sure looping would be the best approach

MIN_DATA1 =data.frame(Strike09 <-  c (1142.757, 1090.971, 1111.644, 1138.833), Strike095 <- c (1206.2435, 1151.5805,1173.4020, 1202.1015), Strike1 <- c (1269.73, 1212.19, 1235.16, 1265.37),Strike105 <- c (1333.216, 1272.800, 1296.918, 1328.639), Strike11 <- c (1396.703,1333.409, 1358.676, 1391.907), Liborrate <- c (0.0505750, 0.0500500, 0.0497078, 0.0496969), TimetoMaturity <- c (0.25, 0.25, 0.25, 0.25), SPX <- c (1269.73, 1212.19, 1235.16, 1265.37), blackscholesc09 <- c(24.995126, 34.905765, 32.103535, 29.686353), blackscholesc095 <- c(37.78425, 50.31239, 45.41761, 43.50957), blackscholesc1 <- c(57.87691, 69.78892, 65.36423, 64.41497), blackscholesc105 <- c(33.47135, 44.10261, 44.12110, 41.11879), blackscholesc11<- c(12.671979, 21.055396, 21.175705, 18.883918))

for (i in 1:length(MIN_DATA)){
results[i] <- optim( x_seed , fn = error_vector_price, gr = NULL,  Strike09[i], Strike095[i], Strike1[i], Strike105[i], Strike11[i], Liborrate[i], TimetoMaturity[h], SPX[i], blackscholesc09[i], blackscholesc095[i], blackscholesc1[i], blackscholesc105[i], blackscholesc11[i], method = "L-BFGS-B", upper = (0.9))
}

results

I hope my problem is clear enough.

Please, provide a reproducible example, so that you maximize your chances of getting some help: [How to make a great R reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). — PaulS, Oct 20 '21 at 15:37
You can do this by nesting, and then using `purrr::map2` or `mapply` (in a straightforward way). We cannot help as you have not provided us with a minimal reproducible example. But look up `tidyr::nest` for this. — Mossa, Oct 20 '21 at 16:14
@Mossa, thank you for you commented. I tried to add a reproducible example as the eventual implementation of your suggestion is not 100% clear to me. Thank you again. — Nazario Maria Cruciano, Oct 20 '21 at 17:29

score 0 · Answer 1 · answered Oct 20 '21 at 18:10

Alright, so the answer is:

results <- vector("list", length = nrow(MIN_DATA1))
with(MIN_DATA1,
     for (i in seq_len(nrow(MIN_DATA1))) {
       results[[i]] <<-
         optim(
           x_seed,
           fn = error_vector_price,
           gr = NULL,
           Strike09 = Strike09[i],
           Strike095 = Strike095[i],
           Strike1 = Strike1[i],
           Strike105 = Strike105[i],
           Strike11 = Strike11[i],
           Liborrate = Liborrate[i],
           TimetoMaturity = TimetoMaturity[i],
           SPX = SPX[i],
           blackscholesc09 = blackscholesc09[i],
           blackscholesc095 = blackscholesc095[i],
           blackscholesc1 = blackscholesc1[i],
           blackscholesc105 = blackscholesc105[i],
           blackscholesc11 = blackscholesc11[i],
           method = "L-BFGS-B",
           upper = (0.9)
         )
     }
)

But there are a bunch of things that could be adjusted with your script. I've made some changes and will paste them here

#' This comes from [SO](https://stackoverflow.com/questions/69648537/how-to-loop-optim-over-a-dataframe-in-r?noredirect=1#comment123111774_69648537)
#'
#'
#'

x_seed <- c(0.64,-0.57,0.25,0.53,0.36)


MIN_DATA1 <-  data.frame(
  Strike09 =
    c(1142.757, 1090.971, 1111.644, 1138.833),
  Strike095 =
    c(1206.2435, 1151.5805, 1173.4020, 1202.1015),
  Strike1 =
    c(1269.73, 1212.19, 1235.16, 1265.37),
  Strike105 =
    c(1333.216, 1272.800, 1296.918, 1328.639),
  Strike11 =
    c(1396.703, 1333.409, 1358.676, 1391.907),
  Liborrate =
    c(0.0505750, 0.0500500, 0.0497078, 0.0496969),
  TimetoMaturity =
    c(0.25, 0.25, 0.25, 0.25),
  SPX =
    c(1269.73, 1212.19, 1235.16, 1265.37),
  blackscholesc09 =
    c(24.995126, 34.905765, 32.103535, 29.686353),
  blackscholesc095 =
    c(37.78425, 50.31239, 45.41761, 43.50957),
  blackscholesc1 =
    c(57.87691, 69.78892, 65.36423, 64.41497),
  blackscholesc105 =
    c(33.47135, 44.10261, 44.12110, 41.11879),
  blackscholesc11 = c(12.671979, 21.055396, 21.175705, 18.883918)
)




MixtureCall <- function(Strike, Liborrate, par_x, TimetoMaturity, SPX) {

  alpha1 <- log(SPX) + (par_x[2]-0.5*(par_x[4]^2)*(par_x[4]^2))*TimetoMaturity
  alpha2 <- log(SPX) + (par_x[3]-0.5*(par_x[5]^2)*(par_x[5]^2))*TimetoMaturity
  beta1  <- (par_x[4]^2)*(TimetoMaturity^0.5)
  beta2  <- (par_x[5]^2)*(TimetoMaturity^0.5)
  theta  <- (par_x[1]^2)/(1+par_x[1]^2)

  D1 <- (-log(Strike) + alpha1 + (beta1 ^ 2)) / beta1
  D2 <- D1 - beta1
  D3 <- (-log(Strike) + alpha2 + (beta2 ^ 2)) / beta2
  D4 <- D3 - beta2
  nD1 <- sapply(D1,  function(x) pnorm(x))
  nD2 <- sapply(D2,  function(x) pnorm(x))
  nD3 <- sapply(D3,  function(x) pnorm(x))
  nD4 <- sapply(D4,  function(x) pnorm(x))

  call1 <- exp(alpha1 + beta1 * beta1 / 2) * nD1 - Strike * nD2
  call2 <- exp(alpha2 + beta2 * beta2 / 2) * nD3 - Strike * nD4

  InstRate <- (log(1 + Liborrate * TimetoMaturity)) / TimetoMaturity

  Mixturecall <- exp(-InstRate * TimetoMaturity) * (theta * call1 + (1 - theta) * call2)

  return(Mixturecall)
}

# MixtureCall(Strike105, Liborrate, par_x, TimetoMaturity, SPX)
# MixtureCall(Strike11, Liborrate, par_x, TimetoMaturity, SPX)

# MixtureCall(Strike = MIN_DATA1$Strike105[1],
#             Liborrate = MIN_DATA1$Liborrate[1],
#             par_x = x_seed,
#             TimetoMaturity = MIN_DATA1$TimetoMaturity[1],
#             SPX = MIN_DATA1$SPX[1])
# MixtureCall(Strike = MIN_DATA1$Strike11[1],
#             Liborrate = MIN_DATA1$Liborrate[1],
#             par_x = x_seed,
#             TimetoMaturity = MIN_DATA1$TimetoMaturity[1],
#             SPX = MIN_DATA1$SPX[1])


MixturePut <- function(Strike, Liborrate, par_x, TimetoMaturity, SPX){

  alpha1 = log(SPX) + (par_x[2]-0.5*(par_x[4]^2)*(par_x[4]^2))*TimetoMaturity
  alpha2 = log(SPX) + (par_x[3]-0.5*(par_x[5]^2)*(par_x[5]^2))*TimetoMaturity
  beta1 = (par_x[4]^2)*(TimetoMaturity^0.5)
  beta2 = (par_x[5]^2)*(TimetoMaturity^0.5)
  theta = (par_x[1]^2)/(1+par_x[1]^2)

  D1 = (-log(Strike) + alpha1 + (beta1 ^ 2)) / beta1
  D2 = D1 - beta1
  D3 = (-log(Strike) + alpha2 + (beta2 ^ 2)) / beta2
  D4 = D3 - beta2
  nD1 = sapply(D1, function(x) pnorm(x))
  nD2 = sapply(D2, function(x) pnorm(x))
  nD3 = sapply(D3, function(x) pnorm(x))
  nD4 = sapply(D4, function(x) pnorm(x))

  call1 = exp(alpha1 + beta1 * beta1 / 2) * nD1 - Strike * nD2
  call2 = exp(alpha2 + beta2 * beta2 / 2) * nD3 - Strike * nD4

  InstRate = (log(1 + Liborrate * TimetoMaturity)) / TimetoMaturity

  Mixturecall = exp(-InstRate * TimetoMaturity) * (theta * call1 + (1 - theta) * call2)

  fwdprice = theta * exp(alpha1 + beta1 * beta1 / 2) + (1 - theta) * exp(alpha2 + beta2 * beta2 / 2)

  Mixtureput = Mixturecall + (Strike - fwdprice) / (1 + Liborrate * TimetoMaturity)

  return(Mixtureput)
}


# MixturePut(Strike09, Liborrate, par_x, TimetoMaturity, SPX)
# MixturePut(Strike095, Liborrate, par_x, TimetoMaturity, SPX)
# MixturePut(Strike1, Liborrate, par_x, TimetoMaturity, SPX)
#
# debug(MixturePut)

# MixturePut(Strike = MIN_DATA1$Strike09[1],
#            Liborrate = MIN_DATA1$Liborrate[1],
#            par_x = x_seed,
#            TimetoMaturity = MIN_DATA1[1],
#            SPX = MIN_DATA1$SPX[1])
# MixturePut(Strike = MIN_DATA1$Strike095[1],
#            Liborrate = MIN_DATA1$Liborrate[1],
#            par_x = x_seed,
#            TimetoMaturity = MIN_DATA1[1],
#            SPX = MIN_DATA1$SPX[1])
# MixturePut(Strike = MIN_DATA1$Strike1[1],
#            Liborrate = MIN_DATA1$Liborrate[1],
#            par_x = x_seed,
#            TimetoMaturity = MIN_DATA1[1],
#            SPX = MIN_DATA1$SPX[1])


error_vector_price <- function(Strike09, Strike095,
                               Strike1, Strike105, Strike11,
                               Liborrate, par_x,
                               TimetoMaturity, SPX,
                               blackscholesc09, blackscholesc095,
                               blackscholesc1, blackscholesc105, blackscholesc11){

  model_price_vector09  <- MixturePut(Strike09, Liborrate, par_x, TimetoMaturity, SPX)
  model_price_vector095 <- MixturePut(Strike095, Liborrate, par_x, TimetoMaturity, SPX)
  model_price_vector1   <- MixturePut(Strike1, Liborrate, par_x, TimetoMaturity, SPX)
  model_price_vector105 <- MixtureCall(Strike105, Liborrate, par_x, TimetoMaturity, SPX)
  model_price_vector11  <- MixtureCall(Strike11, Liborrate, par_x, TimetoMaturity, SPX)

  error_vector_price <-
    (((blackscholesc09 - model_price_vector09) ^ 2 + (blackscholesc095 - model_price_vector095) ^
        2 + (blackscholesc1 - model_price_vector1) ^ 2 + (blackscholesc105 - model_price_vector105) ^
        2 + (blackscholesc11 - model_price_vector11) ^ 2
    ))

  return(error_vector_price)
}



results <-
  optim(
    x_seed ,
    fn = error_vector_price,
    gr = NULL,
    Strike09 = 1142.757,
    Strike095 = 1206.2435,
    Strike1 = 1269.73,
    Strike105 = 1333.216,
    Strike11 = 1396.703,
    Liborrate = 0.0505750,
    TimetoMaturity = 0.25,
    SPX = 1269.73,
    blackscholesc09 = 24.995126,
    blackscholesc095 = 37.78425,
    blackscholesc1 = 57.87691,
    blackscholesc105 = 33.47135,
    blackscholesc11 = 12.671979,
    method = "L-BFGS-B",
    upper = (1)
  )

x_star <- results$par
x_star

results <- vector("list", length = nrow(MIN_DATA1))
with(MIN_DATA1,
     for (i in seq_len(nrow(MIN_DATA1))) {
       results[[i]] <<-
         optim(
           x_seed,
           fn = error_vector_price,
           gr = NULL,
           Strike09 = Strike09[i],
           Strike095 = Strike095[i],
           Strike1 = Strike1[i],
           Strike105 = Strike105[i],
           Strike11 = Strike11[i],
           Liborrate = Liborrate[i],
           TimetoMaturity = TimetoMaturity[i],
           SPX = SPX[i],
           blackscholesc09 = blackscholesc09[i],
           blackscholesc095 = blackscholesc095[i],
           blackscholesc1 = blackscholesc1[i],
           blackscholesc105 = blackscholesc105[i],
           blackscholesc11 = blackscholesc11[i],
           method = "L-BFGS-B",
           upper = (0.9)
         )
     }
)

library(tidyverse)
results %>%
  enframe() %>%
  unnest_wider(value) %>%
  unnest_wider(par)

How to loop optim over a dataframe in R

1 Answers1