# Hauck基金公司案例：非线性规划求解
# 目标：通过选择投资组合投资于每个共有基金的比例，以尽可能地接近标准普尔500的回报

## 加载必要的包----
if (!require("nloptr")) install.packages("nloptr")
library(nloptr)

## 定义目标函数----
eval_f <- function(x) {
  # 提取决策变量
  FS <- x[1]; IB <- x[2]; LG <- x[3]; LV <- x[4]; SG <- x[5]; SV <- x[6]
  
  # 计算每年的投资组合回报
  R1 <- 0.1006*FS + 0.1764*IB + 0.3241*LG + 0.3236*LV + 0.3344*SG + 0.2456*SV
  R2 <- 0.1312*FS + 0.0325*IB + 0.1871*LG + 0.2061*LV + 0.1940*SG + 0.2532*SV
  R3 <- 0.1347*FS + 0.0751*IB + 0.3328*LG + 0.1293*LV + 0.0385*SG - 0.0670*SV
  R4 <- 0.4542*FS - 0.0133*IB + 0.4146*LG + 0.0706*LV + 0.5868*SG + 0.0543*SV
  R5 <- -0.2193*FS + 0.0736*IB - 0.2326*LG - 0.0537*LV - 0.0902*SG + 0.1731*SV
  
  # 计算目标函数值（与标准普尔500收益的差距平方和）
  obj <- (R1 - 0.25)^2 + (R2 - 0.20)^2 + (R3 - 0.08)^2 + 
         (R4 - 0.30)^2 + (R5 + 0.10)^2
  
  return(obj)
}

## 定义目标函数的梯度----
eval_grad_f <- function(x) {
  # 提取决策变量
  FS <- x[1]; IB <- x[2]; LG <- x[3]; LV <- x[4]; SG <- x[5]; SV <- x[6]
  
  # 计算每年的投资组合回报与目标的差距
  R1 <- 0.1006*FS + 0.1764*IB + 0.3241*LG + 0.3236*LV + 0.3344*SG + 0.2456*SV - 0.25
  R2 <- 0.1312*FS + 0.0325*IB + 0.1871*LG + 0.2061*LV + 0.1940*SG + 0.2532*SV - 0.20
  R3 <- 0.1347*FS + 0.0751*IB + 0.3328*LG + 0.1293*LV + 0.0385*SG - 0.0670*SV - 0.08
  R4 <- 0.4542*FS - 0.0133*IB + 0.4146*LG + 0.0706*LV + 0.5868*SG + 0.0543*SV - 0.30
  R5 <- -0.2193*FS + 0.0736*IB - 0.2326*LG - 0.0537*LV - 0.0902*SG + 0.1731*SV + 0.10
  
  # 计算梯度
  grad <- c(
    2*(0.1006*R1 + 0.1312*R2 + 0.1347*R3 + 0.4542*R4 - 0.2193*R5),  # dF/dFS
    2*(0.1764*R1 + 0.0325*R2 + 0.0751*R3 - 0.0133*R4 + 0.0736*R5),  # dF/dIB
    2*(0.3241*R1 + 0.1871*R2 + 0.3328*R3 + 0.4146*R4 - 0.2326*R5),  # dF/dLG
    2*(0.3236*R1 + 0.2061*R2 + 0.1293*R3 + 0.0706*R4 - 0.0537*R5),  # dF/dLV
    2*(0.3344*R1 + 0.1940*R2 + 0.0385*R3 + 0.5868*R4 - 0.0902*R5),  # dF/dSG
    2*(0.2456*R1 + 0.2532*R2 - 0.0670*R3 + 0.0543*R4 + 0.1731*R5)   # dF/dSV
  )
  
  return(grad)
}

## 定义约束条件----
eval_g <- function(x) {
  # 所有投资比例之和等于1的约束
  return(sum(x) - 1)
}

## 定义约束条件的梯度----
eval_grad_g <- function(x) {
  # 约束条件 sum(x) - 1 = 0 的梯度
  return(rep(1, 6))
}

## 设置初始值和边界条件----
x0 <- rep(1/6, 6)  # 初始均匀分配
lb <- rep(0, 6)    # 下界
ub <- rep(1, 6)    # 上界

## 使用 nloptr 求解----
result <- nloptr(
  x0 = x0,  # 初始值
  eval_f = eval_f,  # 目标函数
  eval_grad_f = eval_grad_f,  # 目标函数的梯度
  lb = lb,  # 下界
  ub = ub,  # 上界
  eval_g_eq = eval_g,  # 约束条件
  eval_jac_g_eq = eval_grad_g,  # 约束条件的梯度
  opts = list(
    "algorithm" = "NLOPT_LD_SLSQP",  # 算法
    "xtol_rel" = 1.0e-8,  # 相对误差
    "maxeval" = 1000  # 最大迭代次数
  )
)

## 创建最优解数据框----
optimal_solution <- data.frame(
  Fund = c("外国股票(FS)", "中期债券(IB)", "大市值成长(LG)", 
           "大市值价值(LV)", "小市值成长(SG)", "小市值价值(SV)"),
  Weight = result$solution
)

## 计算最优投资组合下的年度回报----
x_opt <- result$solution
FS <- x_opt[1]; IB <- x_opt[2]; LG <- x_opt[3]
LV <- x_opt[4]; SG <- x_opt[5]; SV <- x_opt[6]

## 计算每年的回报率----
R1 <- 0.1006*FS + 0.1764*IB + 0.3241*LG + 0.3236*LV + 0.3344*SG + 0.2456*SV
R2 <- 0.1312*FS + 0.0325*IB + 0.1871*LG + 0.2061*LV + 0.1940*SG + 0.2532*SV
R3 <- 0.1347*FS + 0.0751*IB + 0.3328*LG + 0.1293*LV + 0.0385*SG - 0.0670*SV
R4 <- 0.4542*FS - 0.0133*IB + 0.4146*LG + 0.0706*LV + 0.5868*SG + 0.0543*SV
R5 <- -0.2193*FS + 0.0736*IB - 0.2326*LG - 0.0537*LV - 0.0902*SG + 0.1731*SV

## 创建年度回报对比表----
annual_returns <- data.frame(
  Year = 1:5,
  Portfolio_Return = 100*c(R1, R2, R3, R4, R5),  # 转换为百分比
  SP500_Return = c(25, 20, 8, 30, -10)  # 标准普尔500收益率（百分比）
)

## 计算目标函数的最优值----
optimal_value <- eval_f(result$solution)*100  # 转换为百分比

# 注释掉打印语句
# print("最优投资组合配置：")
# print(optimal_solution)
# print("\n年度回报对比：")
# print(annual_returns)
# print(paste("目标函数的最优值：", scales::percent(optimal_value, accuracy=0.000001)))