R Markdown
# 加载包
library(quantmod)
## 载入需要的程辑包：xts
## 载入需要的程辑包：zoo
##
## 载入程辑包：'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## 载入需要的程辑包：TTR
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(PerformanceAnalytics)
##
## 载入程辑包：'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
library(CVXR)
##
## 载入程辑包：'CVXR'
## The following object is masked from 'package:stats':
##
## power
library(ggplot2)
# 读取数据（假设数据格式为：日期 + 股票收益率列）
data <- read.csv("C:/Users/L/Desktop/new homework/data(1).csv", header = TRUE, stringsAsFactors = FALSE)
rownames(data) <- as.Date(data$close_price,format = "%Y/%m/%d") # 假设第一列是日期
returns <- data[, -1] # 去掉日期列，保留收益率数据
# 检查列名（确认股票代码列）
print(colnames(data))
## [1] "close_price" "X430047.BJ" "X430090.BJ" "X430139.BJ" "X430198.BJ"
## [6] "X430300.BJ" "X430418.BJ" "X430425.BJ" "X430476.BJ" "X430478.BJ"
## [11] "X430489.BJ" "X430510.BJ" "X430556.BJ" "X430564.BJ" "X430685.BJ"
## [16] "X430718.BJ" "X830799.BJ" "X830809.BJ" "X830832.BJ" "X830839.BJ"
## [21] "X830879.BJ"
# 根据学号后四位2275选择股票，有重复数字2，改为2375（根据实际列名调整）
selected_returns <- data[, c("X430090.BJ", "X430139.BJ", "X430425.BJ", "X430300.BJ")]
# 检查数据
head(selected_returns)
## X430090.BJ X430139.BJ X430425.BJ X430300.BJ
## 2022-01-04 10.03 20.26 24.5 8.86
## 2022-01-05 9.60 19.47 24.7 8.88
## 2022-01-06 8.50 19.49 24.7 8.91
## 2022-01-07 8.29 19.35 24.7 8.67
## 2022-01-10 8.76 19.61 24.7 8.70
## 2022-01-11 8.30 20.00 24.6 8.61
colnames(selected_returns) <- c("Stock2_1", "Stock2_2", "Stock7", "Stock5") # 重命名避免重复
#最小风险权重与直方图
# 定义优化问题（最小化方差）
min_risk_weights <- function(returns) {
n <- ncol(returns)
w <- Variable(n)
risk <- quad_form(w, cov(returns))
prob <- Problem(Minimize(risk), constraints = list(w >= 0, sum(w) == 1))
result <- solve(prob)
return(result$getValue(w))
}
weights_min_risk <- min_risk_weights(selected_returns)
rownames(weights_min_risk) <- colnames(selected_returns)
# 打印权重
print("最小风险权重：")
## [1] "最小风险权重："
print(weights_min_risk)
## [,1]
## Stock2_1 2.015575e-23
## Stock2_2 -5.698232e-23
## Stock7 -1.632472e-24
## Stock5 1.000000e+00
# 绘制直方图
ggplot(data.frame(Stock = colnames(selected_returns), Weight = weights_min_risk[,1]),
aes(x = Stock, y = Weight)) +
geom_bar(stat = "identity", fill = "steelblue") +
labs(title = "Minimum Risk Portfolio Weights", y = "Weight") +
theme_minimal()
#最大收益权重与直方图
max_return_weights <- function(returns) {
n <- ncol(returns)
w <- Variable(n)
expected_return <- t(colMeans(returns)) %*% w
prob <- Problem(Maximize(expected_return), constraints = list(w >= 0, sum(w) == 1))
result <- solve(prob)
return(result$getValue(w))
}
weights_max_return <- max_return_weights(selected_returns)
rownames(weights_max_return) <- colnames(selected_returns)
# 打印权重
print("最大收益权重：")
## [1] "最大收益权重："
print(weights_max_return)
## [,1]
## Stock2_1 -2.121073e-28
## Stock2_2 -2.126238e-28
## Stock7 1.000000e+00
## Stock5 -2.122475e-28
# 绘制直方图
ggplot(data.frame(Stock = colnames(selected_returns), Weight = weights_max_return[,1]),
aes(x = Stock, y = Weight)) +
geom_bar(stat = "identity", fill = "darkred") +
labs(title = "Maximum Return Portfolio Weights", y = "Weight") +
theme_minimal()
#最大夏普比率权重与直方图
max_sharpe_weights <- function(returns, rf = 0) {
n <- ncol(returns)
mu <- colMeans(returns) - rf
Sigma <- cov(returns)
# 使用等效凸优化问题
y <- Variable(n)
prob <- Problem(Minimize(quad_form(y, Sigma)),
constraints = list(t(mu) %*% y == 1, y >= 0))
result <- solve(prob, solver = "ECOS")
if(result$status != "optimal") {
stop("优化失败: ", result$status)
}
# 归一化得到最终权重
w <- result$getValue(y)
return(w / sum(w))
}
# 1. 加载必要包
library(CVXR)
library(Matrix)
# 2. 准备数据
selected_returns <- as.matrix(na.omit(selected_returns))
cov_matrix <- Matrix::nearPD(cov(selected_returns))$mat
# 3. 计算最大夏普比率权重
weights_max_sharpe <- max_sharpe_weights(selected_returns)
# 4. 添加行名
rownames(weights_max_sharpe) <- colnames(selected_returns)
# 5. 输出结果
print("最大夏普比率权重：")
## [1] "最大夏普比率权重："
print(weights_max_sharpe)
## [,1]
## Stock2_1 2.257766e-09
## Stock2_2 1.977207e-09
## Stock7 9.417137e-03
## Stock5 9.905829e-01
# 6. 可视化
barplot(weights_max_sharpe[,1],
main = "最大夏普比率组合权重",
xlab = "资产",
ylab = "权重",
col = "darkgreen")
#关于有效前沿
efficient_frontier <- function(returns, n_points = 50) {
cov_mat <- cov(returns)
mean_ret <- colMeans(returns)
n <- ncol(returns)
# 生成目标收益序列
target_returns <- seq(min(mean_ret), max(mean_ret), length.out = n_points)
target_risks <- numeric(n_points)
weights_list <- list()
for (i in 1:n_points) {
w <- Variable(n)
risk <- quad_form(w, cov_mat)
prob <- Problem(Minimize(risk),
constraints = list(t(mean_ret) %*% w == target_returns[i],
w >= 0, sum(w) == 1))
result <- solve(prob)
target_risks[i] <- sqrt(result$value)
weights_list[[i]] <- as.vector(result$getValue(w)) # 转换为向量
}
# 将权重列表转换为矩阵
weights_matrix <- do.call(rbind, weights_list)
colnames(weights_matrix) <- colnames(returns)
# 返回数据框
data.frame(
Return = target_returns,
Risk = target_risks,
Weights = weights_matrix
)
}
# 生成有效前沿数据（假设selected_returns已存在）
frontier_data <- efficient_frontier(selected_returns)
# 提取关键点
min_risk_point <- frontier_data[which.min(frontier_data$Risk), ]
max_return_point <- frontier_data[which.max(frontier_data$Return), ]
max_sharpe_point <- frontier_data[which.max((frontier_data$Return - 0) / frontier_data$Risk), ]
# 绘制有效前沿
ggplot(frontier_data, aes(x = Risk, y = Return)) +
geom_line(color = "blue") +
geom_point(data = min_risk_point, aes(x = Risk, y = Return), color = "red", size = 3) +
geom_point(data = max_return_point, aes(x = Risk, y = Return), color = "green", size = 3) +
geom_point(data = max_sharpe_point, aes(x = Risk, y = Return), color = "purple", size = 3) +
labs(title = "Efficient Frontier", x = "Risk (Standard Deviation)", y = "Expected Return") +
annotate("text", x = min_risk_point$Risk, y = min_risk_point$Return,
label = "Min Risk", vjust = -1) +
annotate("text", x = max_return_point$Risk, y = max_return_point$Return,
label = "Max Return", vjust = -1) +
annotate("text", x = max_sharpe_point$Risk, y = max_sharpe_point$Return,
label = "Max Sharpe", vjust = -1) +
theme_minimal()
