# Markowitz投资组合模型：考虑期望收益和风险的非线性规划
# 目标：在满足最小期望收益要求下，最小化投资组合的方差

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

## 定义目标函数（投资组合方差）----
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
  
  # 计算期望收益
  R_mean <- (R1 + R2 + R3 + R4 + R5) / 5
  
  # 计算方差
  variance <- ((R1 - R_mean)^2 + (R2 - R_mean)^2 + (R3 - R_mean)^2 + 
               (R4 - R_mean)^2 + (R5 - R_mean)^2) / 5
  
  return(variance)
}

## 定义目标函数的梯度----
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
  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
  
  # 计算期望收益
  R_mean <- (R1 + R2 + R3 + R4 + R5) / 5
  
  # 计算梯度
  grad <- c(
    2/5 * (0.1006*(R1-R_mean) + 0.1312*(R2-R_mean) + 0.1347*(R3-R_mean) + 
           0.4542*(R4-R_mean) - 0.2193*(R5-R_mean)),
    2/5 * (0.1764*(R1-R_mean) + 0.0325*(R2-R_mean) + 0.0751*(R3-R_mean) - 
           0.0133*(R4-R_mean) + 0.0736*(R5-R_mean)),
    2/5 * (0.3241*(R1-R_mean) + 0.1871*(R2-R_mean) + 0.3328*(R3-R_mean) + 
           0.4146*(R4-R_mean) - 0.2326*(R5-R_mean)),
    2/5 * (0.3236*(R1-R_mean) + 0.2061*(R2-R_mean) + 0.1293*(R3-R_mean) + 
           0.0706*(R4-R_mean) - 0.0537*(R5-R_mean)),
    2/5 * (0.3344*(R1-R_mean) + 0.1940*(R2-R_mean) + 0.0385*(R3-R_mean) + 
           0.5868*(R4-R_mean) - 0.0902*(R5-R_mean)),
    2/5 * (0.2456*(R1-R_mean) + 0.2532*(R2-R_mean) - 0.0670*(R3-R_mean) + 
           0.0543*(R4-R_mean) + 0.1731*(R5-R_mean))
  )
  
  return(grad)
}

## 定义约束条件----
eval_g <- 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
  
  # 计算期望收益
  R_mean <- (R1 + R2 + R3 + R4 + R5) / 5
  
  # 返回约束条件
  return(c(
    sum(x) - 1,           # 投资比例之和等于1
    R_mean - 0.10         # 期望收益至少为10%
  ))
}

## 定义约束条件的雅可比矩阵----
eval_jac_g <- function(x) {
  # 提取决策变量
  FS <- x[1]; IB <- x[2]; LG <- x[3]; LV <- x[4]; SG <- x[5]; SV <- x[6]
  
  # 计算期望收益对各变量的偏导数
  # 正确计算：考虑所有五年的收益情况
  dR_mean <- c(
    (0.1006 + 0.1312 + 0.1347 + 0.4542 - 0.2193) / 5,  # FS的偏导数
    (0.1764 + 0.0325 + 0.0751 - 0.0133 + 0.0736) / 5,  # IB的偏导数
    (0.3241 + 0.1871 + 0.3328 + 0.4146 - 0.2326) / 5,  # LG的偏导数
    (0.3236 + 0.2061 + 0.1293 + 0.0706 - 0.0537) / 5,  # LV的偏导数
    (0.3344 + 0.1940 + 0.0385 + 0.5868 - 0.0902) / 5,  # SG的偏导数
    (0.2456 + 0.2532 - 0.0670 + 0.0543 + 0.1731) / 5   # SV的偏导数
  )
  
  # 返回雅可比矩阵
  return(rbind(
    rep(1, 6),           # 投资比例之和的梯度
    dR_mean              # 期望收益的梯度
  ))
}

## 设置初始值和边界条件----
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_jac_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

## 计算期望收益和方差----
R_mean <- (R1 + R2 + R3 + R4 + R5) / 5
variance <- 100*((R1 - R_mean)^2 + (R2 - R_mean)^2 + (R3 - R_mean)^2 + 
             (R4 - R_mean)^2 + (R5 - R_mean)^2) / 5

## 创建年度回报表----
annual_returns <- data.frame(
  Year = 1:5,
  Portfolio_Return = c(R1, R2, R3, R4, R5)  # 收益水平
  )|>  
  # 更改year为character
  mutate(Year = as.character(Year)) |>
  # 最后一行添加平均收益
  bind_rows(
    data.frame(Year = "平均收益R", Portfolio_Return = R_mean)
  )

## 注释掉打印语句----
# print("最优投资组合配置（权重）：")
# print(optimal_solution)
# print("\n年度回报：")
# print(annual_returns |> mutate(Portfolio_Return =  scales::percent(Portfolio_Return, accuracy=0.00001)))
# print(paste("期望收益率：", scales::percent(R_mean, accuracy=0.01)))
# print(paste(最优目标函数值（投资组合方差）：", scales::percent(variance, accuracy=0.01))) "