Untitled

# 加载包
library(readxl)
library(quadprog)
library(ggplot2)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

# 读取数据
df <- read_excel("data.xls", col_names = FALSE)

## New names:
## • `` -> `...1`
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`

# 提取股票代码和名称
stock_codes <- as.character(df[3, 2:21])
stock_names <- as.character(df[4, 2:21])

# 学号2934选中股票
selected_indices <- c(2, 9, 3, 4)
selected_codes <- stock_codes[selected_indices]
selected_names <- stock_names[selected_indices]

cat("选中股票:\n")

## 选中股票:

for(i in 1:4) {
  cat(sprintf("%d: %s - %s\n", selected_indices[i], selected_codes[i], selected_names[i]))
}

## 2: 10.029999999999999 - 9.5999999999999996
## 9: 21.989999999999998 - 21.510000000000002
## 3: 20.260000000000002 - 19.469999999999999
## 4: 9.6999999999999993 - 9.7699999999999996

# 提取价格数据
price_data <- df[5:nrow(df), c(1, selected_indices + 1)]
colnames(price_data) <- c("Date", selected_codes)

# 转换格式
price_data$Date <- as.Date(as.numeric(price_data$Date), origin = "1899-12-30")
for(col in selected_codes) {
  price_data[[col]] <- as.numeric(price_data[[col]])
}
price_data <- na.omit(price_data)
# 计算对数收益率
returns <- price_data[, -1]
for(col in names(returns)) {
  returns[[col]] <- c(NA, diff(log(returns[[col]])))
}
returns <- na.omit(returns)

# 年化参数
trading_days <- 252
mean_returns <- colMeans(returns) * trading_days
cov_matrix <- cov(returns) * trading_days

cat("年化收益率:\n")

## 年化收益率:

print(round(mean_returns, 4))

## 10.029999999999999 21.989999999999998 20.260000000000002 9.6999999999999993 
##            -0.9444            -0.1215            -0.5093            -0.5984

# 最小方差组合
n <- length(mean_returns)
Dmat <- 2 * cov_matrix
dvec <- rep(0, n)
Amat <- cbind(rep(1, n), diag(n))
bvec <- c(1, rep(0, n))

min_risk_result <- solve.QP(Dmat, dvec, Amat, bvec, meq = 1)
min_risk_weights <- min_risk_result$solution

min_risk_return <- sum(min_risk_weights * mean_returns)
min_risk_risk <- sqrt(sum(min_risk_weights * (cov_matrix %*% min_risk_weights)))

cat("最小风险组合权重:\n")

## 最小风险组合权重:

for(i in 1:4) {
  cat(sprintf("  %s: %.4f\n", selected_names[i], min_risk_weights[i]))
}

##   9.5999999999999996: 0.0066
##   21.510000000000002: 0.0304
##   19.469999999999999: 0.2659
##   9.7699999999999996: 0.6970

cat(sprintf("\n收益: %.4f (%.2f%%)\n", min_risk_return, min_risk_return*100))

## 
## 收益: -0.5625 (-56.25%)

cat(sprintf("风险: %.4f (%.2f%%)\n", min_risk_risk, min_risk_risk*100))

## 风险: 0.2316 (23.16%)

max_return_idx <- which.max(mean_returns)
max_return_weights <- rep(0, 4)
max_return_weights[max_return_idx] <- 1

max_return_return <- sum(max_return_weights * mean_returns)
max_return_risk <- sqrt(sum(max_return_weights * (cov_matrix %*% max_return_weights)))

cat("最大收益组合权重:\n")

## 最大收益组合权重:

for(i in 1:4) {
  cat(sprintf("  %s: %.4f\n", selected_names[i], max_return_weights[i]))
}

##   9.5999999999999996: 0.0000
##   21.510000000000002: 1.0000
##   19.469999999999999: 0.0000
##   9.7699999999999996: 0.0000

cat(sprintf("\n收益: %.4f (%.2f%%)\n", max_return_return, max_return_return*100))

## 
## 收益: -0.1215 (-12.15%)

cat(sprintf("风险: %.4f (%.2f%%)\n", max_return_risk, max_return_risk*100))

## 风险: 1.1565 (115.65%)

cat(sprintf("选中股票: %s\n", selected_names[max_return_idx]))

## 选中股票: 21.510000000000002

# 生成有效前沿
generate_efficient_frontier <- function(mean_returns, cov_matrix, n_points = 100) {
  n <- length(mean_returns)
  target_returns <- seq(min(mean_returns), max(mean_returns), length.out = n_points)
  
  ef_risks <- c()
  ef_returns <- c()
  ef_weights <- list()
  
  for(target in target_returns) {
    Dmat <- 2 * cov_matrix
    dvec <- rep(0, n)
    Amat <- cbind(rep(1, n), mean_returns, diag(n))
    bvec <- c(1, target, rep(0, n))
    
    result <- tryCatch({
      solve.QP(Dmat, dvec, Amat, bvec, meq = 2)
    }, error = function(e) NULL)
    
    if(!is.null(result)) {
      w <- result$solution
      risk <- sqrt(sum(w * (cov_matrix %*% w)))
      ef_risks <- c(ef_risks, risk)
      ef_returns <- c(ef_returns, target)
      ef_weights[[length(ef_weights) + 1]] <- w
    }
  }
  
  return(list(risk = ef_risks, return = ef_returns, weights = ef_weights))
}

ef_data <- generate_efficient_frontier(mean_returns, cov_matrix)
sharpe_ratios <- ef_data$return / ef_data$risk
max_sharpe_idx <- which.max(sharpe_ratios)

max_sharpe_weights <- ef_data$weights[[max_sharpe_idx]]
max_sharpe_return <- ef_data$return[max_sharpe_idx]
max_sharpe_risk <- ef_data$risk[max_sharpe_idx]

cat("最大夏普比率组合权重:\n")

## 最大夏普比率组合权重:

for(i in 1:4) {
  cat(sprintf("  %s: %.4f\n", selected_names[i], max_sharpe_weights[i]))
}

##   9.5999999999999996: 0.0000
##   21.510000000000002: 1.0000
##   19.469999999999999: -0.0000
##   9.7699999999999996: 0.0000

cat(sprintf("\n收益: %.4f (%.2f%%)\n", max_sharpe_return, max_sharpe_return*100))

## 
## 收益: -0.1215 (-12.15%)

cat(sprintf("风险: %.4f (%.2f%%)\n", max_sharpe_risk, max_sharpe_risk*100))

## 风险: 1.1565 (115.65%)

cat(sprintf("夏普比率: %.4f\n", sharpe_ratios[max_sharpe_idx]))

## 夏普比率: -0.1051

# 生成随机投资组合
set.seed(123)
n_random <- 3000
random_risks <- numeric(n_random)
random_returns <- numeric(n_random)
random_sharpe <- numeric(n_random)

for(i in 1:n_random) {
  w <- runif(4)
  w <- w / sum(w)
  random_returns[i] <- sum(w * mean_returns)
  random_risks[i] <- sqrt(sum(w * (cov_matrix %*% w)))
  random_sharpe[i] <- random_returns[i] / random_risks[i]
}

# 创建数据框
random_df <- data.frame(risk = random_risks, return = random_returns, sharpe = random_sharpe)

# 特殊点
special_points <- data.frame(
  name = c("Min Risk", "Max Return", "Max Sharpe"),
  risk = c(min_risk_risk, max_return_risk, max_sharpe_risk),
  return = c(min_risk_return, max_return_return, max_sharpe_return)
)

# 绘图
p <- ggplot() +
  geom_point(data = random_df, aes(x = risk, y = return, color = sharpe), 
             alpha = 0.5, size = 1.5) +
  scale_color_gradient(low = "purple", high = "yellow", name = "Sharpe") +
  geom_line(data = as.data.frame(ef_data), aes(x = risk, y = return), 
            color = "red", linewidth = 1.5) +
  geom_point(data = special_points, aes(x = risk, y = return), 
             size = 5, shape = 21, fill = "white", stroke = 2, color = "black") +
  geom_text(data = special_points, aes(x = risk, y = return, label = name),
            vjust = -1.5, size = 4, fontface = "bold", color = "black") +
  labs(title = "Portfolio Optimization & Efficient Frontier",
       subtitle = "Student ID: 2934 | X: Risk, Y: Return",
       x = "Risk (Standard Deviation)", 
       y = "Expected Return") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 11, hjust = 0.5),
    axis.title = element_text(size = 12),
    legend.position = "right",
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color = "black", fill = NA, linewidth = 0.5)
  )

print(p)

results_summary <- data.frame(
  组合类型 = c("最小风险", "最大收益", "最大夏普"),
  同辉信息 = round(c(min_risk_weights[1], max_return_weights[1], max_sharpe_weights[1]), 4),
  峆一药业 = round(c(min_risk_weights[2], max_return_weights[2], max_sharpe_weights[2]), 4),
  华岭股份 = round(c(min_risk_weights[3], max_return_weights[3], max_sharpe_weights[3]), 4),
  微创光电 = round(c(min_risk_weights[4], max_return_weights[4], max_sharpe_weights[4]), 4),
  收益率 = round(c(min_risk_return, max_return_return, max_sharpe_return) * 100, 2),
  风险 = round(c(min_risk_risk, max_return_risk, max_sharpe_risk) * 100, 2),
  夏普比率 = round(c(min_risk_return/min_risk_risk, max_return_return/max_return_risk, sharpe_ratios[max_sharpe_idx]), 4)
)

print(results_summary, row.names = FALSE)

##  组合类型 同辉信息 峆一药业 华岭股份 微创光电 收益率   风险 夏普比率
##  最小风险   0.0066   0.0304   0.2659    0.697 -56.25  23.16  -2.4281
##  最大收益   0.0000   1.0000   0.0000    0.000 -12.15 115.65  -0.1051
##  最大夏普   0.0000   1.0000   0.0000    0.000 -12.15 115.65  -0.1051

# 1. 最小风险组合权重直方图
p_min_risk <- data.frame(
  股票 = selected_names,
  权重 = min_risk_weights
) %>%
  ggplot(aes(x = 股票, y = 权重, fill = 股票)) +
  geom_col(alpha = 0.8, color = "black") +
  geom_text(aes(label = round(权重, 4)), vjust = -0.3, size = 4, fontface = "bold") +
  scale_fill_brewer(palette = "Set2") +
  labs(title = "最小风险组合 - 资产权重分布",
       x = "股票名称", y = "投资权重") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        legend.position = "none")

print(p_min_risk)

# 2. 最大收益组合权重直方图
p_max_return <- data.frame(
  股票 = selected_names,
  权重 = max_return_weights
) %>%
  ggplot(aes(x = 股票, y = 权重, fill = 股票)) +
  geom_col(alpha = 0.8, color = "black") +
  geom_text(aes(label = round(权重, 4)), vjust = -0.3, size = 4, fontface = "bold") +
  scale_fill_brewer(palette = "Set1") +
  labs(title = "最大收益组合 - 资产权重分布",
       x = "股票名称", y = "投资权重") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        legend.position = "none")

print(p_max_return)

# 3. 最大夏普比率组合权重直方图
p_max_sharpe <- data.frame(
  股票 = selected_names,
  权重 = max_sharpe_weights
) %>%
  ggplot(aes(x = 股票, y = 权重, fill = 股票)) +
  geom_col(alpha = 0.8, color = "black") +
  geom_text(aes(label = round(权重, 4)), vjust = -0.3, size = 4, fontface = "bold") +
  scale_fill_brewer(palette = "Set3") +
  labs(title = "最大夏普比率组合 - 资产权重分布",
       x = "股票名称", y = "投资权重") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
        legend.position = "none")

print(p_max_sharpe)