基于长面板数据的R语言投资组合优化问题求助

核心错误原因 #

1. 变量名不一致：代码混用return_dat和return_data，导致筛选资产时无数据返回，最终asset_returns列数为0，触发维度不兼容错误。
2. 行业标签提取错误：原数据为行业虚拟变量（Sector1/Sector2），但代码直接调用不存在的Sector列，导致行业列表为空。
3. 无效重复抽样：sample(..., replace=TRUE)会重复选中同一行业/资产，既违背跨行业分散需求，也可能导致可用资产数量不足。
4. 奇异矩阵风险：当选中资产完全相关时，协方差矩阵奇异，会导致二次规划求解失败。

修正后的完整代码 #

library(tidyverse)
library(xts)
library(quadprog)

# ---------------------- 1. 数据预处理：生成行业标签 ----------------------
# 从虚拟变量提取每个公司对应的行业（根据实际行业数量扩展case_when）
return_dat <- return_dat %>%
  mutate(Sector = case_when(
    Sector1 == 1 ~ "Sector1",
    Sector2 == 1 ~ "Sector2",
    # 添加其他行业规则
    TRUE ~ NA_character_
  )) %>%
  filter(!is.na(Sector))  # 移除无行业标签的公司

# ---------------------- 2. 定义参数 ----------------------
assets <- unique(return_dat$CompanyID)
asset_min <- 5
asset_max <- 205
asset_step <- 5
asset_range <- seq(from = asset_min, to = asset_max, by = asset_step)

sectors <- unique(return_dat$Sector)
sector_min <- 2
sector_max <- min(10, length(sectors))  # 最大行业数不超过实际存在的行业数
sector_step <- 2
sector_range <- seq(from = sector_min, to = sector_max, by = sector_step)

risk_free_rate <- 0.02  # 定义无风险利率，可根据实际调整

# 初始化结果数据框（指定列类型，避免类型混乱）
Results <- data.frame(
  num_sectors = integer(),
  num_assets = integer(),
  return = numeric(),
  risk = numeric(),
  sharpe_ratio = numeric()
)

# ---------------------- 3. 循环优化 ----------------------
for (k in sector_range) {
  # 随机抽取k个不同行业
  sector_sample <- sample(sectors, k, replace = FALSE)
  asset_dat <- return_dat %>% filter(Sector %in% sector_sample)
  
  # 检查该行业组合下的可用资产数量
  available_assets <- unique(asset_dat$CompanyID)
  if (length(available_assets) < asset_min) next  # 可用资产不足，跳过
  
  for (i in asset_range) {
    # 跳过超过可用资产数量的情况
    if (i > length(available_assets)) next
    
    # 随机抽取i个不同资产
    asset_sample <- sample(available_assets, i, replace = FALSE)
    asset_subset <- asset_dat %>% filter(CompanyID %in% asset_sample)
    
    # 转换为宽格式（资产列，日期行）
    asset_returns_wide <- asset_subset %>%
      select(CompanyID, Date, PR) %>%
      pivot_wider(names_from = CompanyID, values_from = PR)
    
    # 转换为xts对象（可选，也可直接用data.frame计算）
    asset_returns <- xts(asset_returns_wide[,-1], order.by = as.Date(asset_returns_wide$Date))
    
    # 处理缺失值：仅保留所有资产都有数据的日期（或根据需求用na.locf填充）
    asset_returns_clean <- na.omit(asset_returns)
    if (nrow(asset_returns_clean) < 12) next  # 数据量太少，跳过（至少12个月数据）
    
    # 计算年化预期收益和协方差矩阵
    exp_returns <- colMeans(asset_returns_clean) * 12
    cov_matrix <- cov(asset_returns_clean, use = "pairwise.complete.obs") * 12
    
    # 检查协方差矩阵是否奇异（避免solve.QP报错）
    if (det(cov_matrix) < 1e-10) next
    
    # ---------------------- 4. 二次规划设置（最小化方差，即最大分散化） ----------------------
    n_assets <- ncol(asset_returns_clean)
    Dmat <- 2 * cov_matrix  # 目标函数是(1/2)w'Cov w，所以乘以2
    dvec <- rep(0, n_assets)  # 最小化方差时dvec为0
    
    # 约束条件：
    # 1. 权重和为1（等式约束）
    # 2. 权重非负（不等式约束）
    Amat <- cbind(rep(1, n_assets), diag(n_assets))  # n_assets行，n_assets+1列
    bvec <- c(1, rep(0, n_assets))
    meq <- 1  # 第一个约束是等式
    
    # 求解二次规划
    solution <- solve.QP(Dmat, dvec, Amat, bvec, meq)
    weights <- solution$solution
    
    # 计算组合指标
    portfolio_return <- sum(exp_returns * weights)
    portfolio_risk <- sqrt(t(weights) %*% cov_matrix %*% weights)
    sharpe_ratio <- (portfolio_return - risk_free_rate) / portfolio_risk
    
    # 保存结果
    Results <- rbind(Results, data.frame(
      num_sectors = k,
      num_assets = i,
      return = portfolio_return,
      risk = portfolio_risk,
      sharpe_ratio = sharpe_ratio
    ))
  }
}

# 查看结果
print(Results)

# 可选：可视化不同行业/资产数量下的夏普比率
ggplot(Results, aes(x = num_assets, y = sharpe_ratio, color = factor(num_sectors))) +
  geom_line() +
  labs(x = "资产数量", y = "夏普比率", color = "行业数量") +
  theme_minimal()

额外说明：最大化夏普比率的调整 #

如果需要最大化夏普比率（而非最小化方差），可使用PortfolioAnalytics包简化优化逻辑：

library(PortfolioAnalytics)

# 示例：针对某一资产子集最大化夏普比率
port_spec <- portfolio.spec(assets = colnames(asset_returns_clean))
port_spec <- add.constraint(port_spec, type = "full_investment")
port_spec <- add.constraint(port_spec, type = "long_only")
port_spec <- add.objective(port_spec, type = "risk", name = "StdDev")
port_spec <- add.objective(port_spec, type = "return", name = "mean", multiplier = 1)

# 求解最大化夏普比率的组合
opt_sharpe <- optimize.portfolio(asset_returns_clean, port_spec, 
                                 optimize_method = "ROI",
                                 trace = FALSE)
weights_sharpe <- extractWeights(opt_sharpe)

内容的提问来源于stack exchange，提问作者Moataz