R语言如何对dataframe每行循环执行optim优化求解最优参数

原代码问题点 #

• 循环范围错误：1:length(MIN_DATA)返回的是数据集列数，应该改为取行数1:nrow(MIN_DATA1)
• 笔误问题：TimetoMaturity[h]中的h未定义，应为当前行索引i
• 参数匹配错误：optim传参未显式指定名称，易出现位置匹配错误
• 上限参数不匹配：L-BFGS-B方法要求上下限长度和待优化参数数量一致，原代码仅传入1个上限值，5个待优化参数需要传入长度为5的上限向量
• 结果存储错误：optim返回的是列表，不能直接用向量/普通数组存储

可行实现方案 #

首先修正数据框定义（原定义中使用<-会导致列名异常，改为=）：

# 目标函数保持你提供的原有代码不变
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)

# 修正数据集定义
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)
)

方案1：修正后循环实现 #

# 初始化结果表，存储5个最优参数+最小误差
results_df <- data.frame(matrix(nrow = nrow(MIN_DATA1), ncol = 6))
colnames(results_df) <- c("par1","par2","par3","par4","par5","min_error")

for (i in 1:nrow(MIN_DATA1)){
  row_data <- MIN_DATA1[i,]
  opt_res <- optim(
    par = x_seed, 
    fn = error_vector_price, 
    Strike09 = row_data$Strike09,
    Strike095 = row_data$Strike095,
    Strike1 = row_data$Strike1,
    Strike105 = row_data$Strike105,
    Strike11 = row_data$Strike11,
    Liborrate = row_data$Liborrate,
    TimetoMaturity = row_data$TimetoMaturity,
    SPX = row_data$SPX,
    blackscholesc09 = row_data$blackscholesc09,
    blackscholesc095 = row_data$blackscholesc095,
    blackscholesc1 = row_data$blackscholesc1,
    blackscholesc105 = row_data$blackscholesc105,
    blackscholesc11 = row_data$blackscholesc11,
    method = "L-BFGS-B", 
    upper = rep(0.9,5) # 5个参数分别设置上限
  )
  results_df[i,1:5] <- opt_res$par
  results_df[i,6] <- opt_res$value
}

# 查看结果
print(results_df)

方案2：purrr向量化实现（更简洁，无需手动写循环） #

library(purrr)
results_df <- pmap_dfr(MIN_DATA1, function(...){
  current_args <- list(...)
  opt_res <- optim(
    par = x_seed,
    fn = error_vector_price,
    Strike09 = current_args$Strike09,
    Strike095 = current_args$Strike095,
    Strike1 = current_args$Strike1,
    Strike105 = current_args$Strike105,
    Strike11 = current_args$Strike11,
    Liborrate = current_args$Liborrate,
    TimetoMaturity = current_args$TimetoMaturity,
    SPX = current_args$SPX,
    blackscholesc09 = current_args$blackscholesc09,
    blackscholesc095 = current_args$blackscholesc095,
    blackscholesc1 = current_args$blackscholesc1,
    blackscholesc105 = current_args$blackscholesc105,
    blackscholesc11 = current_args$blackscholesc11,
    method = "L-BFGS-B",
    upper = rep(0.9,5)
  )
  return(data.frame(
    par1 = opt_res$par[1],
    par2 = opt_res$par[2],
    par3 = opt_res$par[3],
    par4 = opt_res$par[4],
    par5 = opt_res$par[5],
    min_error = opt_res$value
  ))
})

# 查看结果
print(results_df)

注意事项 #

如果不同参数的上限要求不同，可以自行修改upper参数的向量值，不需要所有参数都用相同上限
如果出现优化不收敛的情况，可以调整初始值x_seed，或者更换优化方法

内容的提问来源于stack exchange，提问作者Nazario Maria Cruciano