一、引言

本次作业基于给定的股票历史数据,运用投资组合理论对资产配置进行分析。根据学号后四位(4437),选取了对应的四只股票(华岭股份、微创光电、辰光医疗、乐创技术)进行研究。
  作业要求:
  1. 求最小风险权重并画出对应直方图
  2. 求最大收益组合的权重并画出直方图
  3. 求夏普比率最大的组合权重并画出直方图
  4. 做投资组合优化并画出有效前沿
  根据以上要求,我们采用了以下工具包:

library(PerformanceAnalytics)
library(quantmod)
library(xts)
library(dygraphs)
library(tibble)
library(tidyr)
library(ggplot2)
library(plotly)
library(timetk)
library(forcats)
library(readxl)
library(dplyr)
library(scales)

二、数据说明

本次使用的数据来自老师给我们的股票历史价格数据,其中: 第4行为变量名称, 第1列为时间(格式为 yyyy/mm/dd), 其余列为各股票价格。
  在数据处理过程中,主要进行了数据读取与清洗、转换为时间序列数据(xts格式)、计算对数收益率、去除缺失值的工作,以下为数据处理代码:

# 读取数据,跳过前4行,第5行开始是数据
data_raw <- read_excel("E:/wangwenyang/data.xls", sheet = "Sheet1", skip = 4, col_names = FALSE)
# 查看数据结构
cat("数据维度:", dim(data_raw), "\n")
# 处理日期(第一列)
dates_char <- as.character(data_raw[[1]])
dates <- as.Date(dates_char, format = "%Y-%m-%d")
dates_posix <- as.POSIXct(dates)
cat("日期范围:", format(min(dates), "%Y-%m-%d"), "至", format(max(dates), "%Y-%m-%d"), "\n")
cat("观测天数:", length(dates), "\n")
# 提取价格数据(第4,5,6,8列)
price_cols <- data_raw[, c(4, 5, 6, 8)]
price_matrix <- as.matrix(price_cols)
mode(price_matrix) <- "numeric"
# 创建xts对象,使用英文列名
price <- xts(price_matrix, order.by = dates_posix)
colnames(price) <- c("Hualing", "Weichuang", "Chenguang", "Lechuang")
# 移除缺失值
price <- na.omit(price)
cat("最终数据范围:", format(start(price)), "至", format(end(price)), "\n")
cat("最终观测数量:", nrow(price), "\n")
# 查看前几行(为了验证数据读取是否正确,输出数据基本信息如下:)
head(price)

2.1 股票价格走势

下图展示了四只股票在样本期间内的价格变化情况,可以观察到不同股票的价格整体上呈现下降态势,而具体的波动特征存在明显差异。

# 创建中文显示的数据用于绘图
price_cn <- price
colnames(price_cn) <- c("华岭股份", "微创光电", "辰光医疗", "乐创技术")

dygraph(price_cn) %>%
  dyOptions(colors = c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728")) %>%
  dyLegend(show = "always", width = 400) %>%
  dyRangeSelector() %>%
  dyAxis("y", label = "价格") %>%
  dyAxis("x", label = "日期")

2.2 收益率走势

为了更准确地衡量资产表现,本文采用对数收益率进行分析。收益率序列能够更好地反映资产的波动性特征。

# 计算对数日收益率
rt_daily <- na.omit(Return.calculate(price, method = "log"))
# 收益率走势图
rt_cn <- rt_daily
colnames(rt_cn) <- c("华岭股份", "微创光电", "辰光医疗", "乐创技术")
dygraph(rt_cn) %>%
  dyOptions(colors = c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728")) %>%
  dyLegend(show = "always", width = 400) %>%
  dyRangeSelector() %>%
  dyAxis("y", label = "对数收益率")

图中可见四只股票的对数收益率均表现出显著的波动聚集性,其中乐创技术与华岭股份的极端波动幅度更大,乐创技术更是出现大幅正负双向跳升,体现出较高的收益风险;而辰光医疗波动相对更为温和,四只股票的差异化波动特征进一步印证了各标的自身的风险收益属性差异。

三、统计量计算

为了更精确地刻画各股票的收益与风险特征,本文进一步计算其平均收益率、协方差矩阵以及标准差等统计指标。

# 日平均收益率
mean_ret_daily <- colMeans(rt_daily)
cat("\n=== 日平均收益率 ===\n")
## 
## === 日平均收益率 ===
print(round(mean_ret_daily, 6))
##   Hualing Weichuang Chenguang  Lechuang 
## -0.002133 -0.002373 -0.001261 -0.002139
# 年化平均收益率 (252个交易日)
mean_ret_annual <- mean_ret_daily * 252
cat("\n=== 年化平均收益率 ===\n")
## 
## === 年化平均收益率 ===
names(mean_ret_annual) <- c("华岭股份", "微创光电", "辰光医疗", "乐创技术")
print(round(mean_ret_annual, 4))
## 华岭股份 微创光电 辰光医疗 乐创技术 
##  -0.5376  -0.5979  -0.3179  -0.5391
# 协方差矩阵 (年化)
cov_mat <- cov(rt_daily) * 252
cat("\n=== 年化协方差矩阵 ===\n")
## 
## === 年化协方差矩阵 ===
colnames(cov_mat) <- c("华岭", "微创", "辰光", "乐创")
rownames(cov_mat) <- c("华岭", "微创", "辰光", "乐创")
print(round(cov_mat, 4))
##        华岭    微创   辰光    乐创
## 华岭 0.1414  0.0222 0.0105  0.0154
## 微创 0.0222  0.0667 0.0228 -0.0150
## 辰光 0.0105  0.0228 0.0950  0.0018
## 乐创 0.0154 -0.0150 0.0018  0.5691
# 年化标准差(风险)
std_annual <- sqrt(diag(cov_mat))
names(std_annual) <- c("华岭股份", "微创光电", "辰光医疗", "乐创技术")
cat("\n=== 年化标准差(风险)===\n")
## 
## === 年化标准差(风险)===
print(round(std_annual, 4))
## 华岭股份 微创光电 辰光医疗 乐创技术 
##   0.3760   0.2582   0.3083   0.7544

四、蒙特卡洛模拟

为了寻找最优投资组合,本文采用蒙特卡洛模拟方法,通过随机生成大量投资组合权重,并计算其对应的收益与风险,从而构建有效前沿。

4.1 模拟设置

在模拟过程中,设定随机生成10000个投资组合,并假设无风险利率为0。通过随机生成权重并进行归一化处理,使得每个组合的权重之和为1,同时,为了保证结果的可重复性,设置随机种子。

4.2 执行模拟

set.seed(4437)  # 设置随机种子

tick <- c("Hualing", "Weichuang", "Chenguang", "Lechuang")
num_port <- 10000

# 存储矩阵和向量
all_wts <- matrix(nrow = num_port, ncol = 4)
port_returns <- vector('numeric', length = num_port)
port_risk <- vector('numeric', length = num_port)
sharpe_ratio <- vector('numeric', length = num_port)

cat("\n开始蒙特卡洛模拟...\n")
## 
## 开始蒙特卡洛模拟...
for (i in 1:num_port) {
  # 随机生成权重
  wts <- runif(4)
  wts <- wts / sum(wts)
  all_wts[i, ] <- wts
  
  # 组合收益(年化)
  port_ret <- sum(wts * mean_ret_daily)
  port_ret <- ((port_ret + 1)^252) - 1
  port_returns[i] <- port_ret
  
  # 组合风险
  port_sd <- sqrt(t(wts) %*% (cov_mat %*% wts))
  port_risk[i] <- port_sd
  
  # 夏普比率(假设无风险利率=0)
  sharpe_ratio[i] <- port_ret / port_sd
  
  # 显示进度
  if(i %% 1000 == 0) cat("已完成", i, "个组合\n")
}
## 已完成 1000 个组合
## 已完成 2000 个组合
## 已完成 3000 个组合
## 已完成 4000 个组合
## 已完成 5000 个组合
## 已完成 6000 个组合
## 已完成 7000 个组合
## 已完成 8000 个组合
## 已完成 9000 个组合
## 已完成 10000 个组合
# 创建结果数据表
portfolio_values <- tibble(
  Hualing = all_wts[,1],
  Weichuang = all_wts[,2],
  Chenguang = all_wts[,3],
  Lechuang = all_wts[,4],
  Return = port_returns,
  Risk = port_risk,
  SharpeRatio = sharpe_ratio
)

cat("\n模拟完成!生成了", nrow(portfolio_values), "个随机组合\n")
## 
## 模拟完成!生成了 10000 个随机组合

五、最优投资组合分析

在完成蒙特卡洛模拟后,可以从所有随机生成的投资组合中筛选出不同优化目标下的最优组合,包括最小风险组合、最大收益组合以及最大夏普比率组合。不同组合代表了不同的投资偏好。

5.1 最小风险组合

最小风险组合是指在所有投资组合中具有最低波动率(标准差)的组合,该组合强调风险最小化,适用于风险厌恶型投资者。

# 找出最小风险组合
min_var <- portfolio_values[which.min(portfolio_values$Risk), ]

cat("\n=== 最小风险组合 ===\n")
## 
## === 最小风险组合 ===
cat("年化风险:", round(min_var$Risk, 4), "\n")
## 年化风险: 0.1996
cat("年化收益:", round(min_var$Return, 4), "\n")
## 年化收益: -0.393
cat("夏普比率:", round(min_var$SharpeRatio, 4), "\n")
## 夏普比率: -1.9691
cat("\n权重分配:\n")
## 
## 权重分配:
min_weights <- round(c(min_var$Hualing, min_var$Weichuang, min_var$Chenguang, min_var$Lechuang), 4)
names(min_weights) <- c("华岭股份", "微创光电", "辰光医疗", "乐创技术")
print(min_weights)
## 华岭股份 微创光电 辰光医疗 乐创技术 
##   0.1727   0.4505   0.3012   0.0757
# 最小风险组合权重直方图
p1 <- data.frame(
  Asset = c("华岭股份", "微创光电", "辰光医疗", "乐创技术"),
  Weights = min_weights
) %>%
  mutate(Asset = factor(Asset, levels = Asset[order(Weights)])) %>%
  ggplot(aes(x = Asset, y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = '资产', y = '权重', title = "最小风险组合权重") +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  geom_text(aes(label = paste0(round(Weights * 100, 1), "%")), 
            vjust = -0.5, size = 4) +
  theme(legend.position = "bottom")

ggplotly(p1) %>% layout(title = "最小风险组合权重分配")

最终我们发现最小风险组合的权重分配是:华岭股份 0.1727 微创光电 0.4505 辰光医疗 0.3012 乐创技术 0.0757。

5.2 最大收益组合

最大收益组合是指在所有投资组合中年化收益率最高的组合,强调收益最大化。

# 找出最大收益组合
max_return <- portfolio_values[which.max(portfolio_values$Return), ]

cat("\n=== 最大收益组合 ===\n")
## 
## === 最大收益组合 ===
cat("年化收益:", round(max_return$Return, 4), "\n")
## 年化收益: -0.2903
cat("年化风险:", round(max_return$Risk, 4), "\n")
## 年化风险: 0.2807
cat("夏普比率:", round(max_return$SharpeRatio, 4), "\n")
## 夏普比率: -1.0342
cat("\n权重分配:\n")
## 
## 权重分配:
max_weights <- round(c(max_return$Hualing, max_return$Weichuang, max_return$Chenguang, max_return$Lechuang), 4)
names(max_weights) <- c("华岭股份", "微创光电", "辰光医疗", "乐创技术")
print(max_weights)
## 华岭股份 微创光电 辰光医疗 乐创技术 
##   0.0463   0.0108   0.8906   0.0523
# 最大收益组合权重直方图
p2 <- data.frame(
  Asset = c("华岭股份", "微创光电", "辰光医疗", "乐创技术"),
  Weights = max_weights
) %>%
  mutate(Asset = factor(Asset, levels = Asset[order(Weights)])) %>%
  ggplot(aes(x = Asset, y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = '资产', y = '权重', title = "最大收益组合权重") +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  geom_text(aes(label = paste0(round(Weights * 100, 1), "%")), 
            vjust = -0.5, size = 4) +
  theme(legend.position = "bottom")

ggplotly(p2) %>% layout(title = "最大收益组合权重分配")

最终我们发现最大收益组合的权重分配是:华岭股份 0.0463 微创光电 0.0108 辰光医疗 0.8906 乐创技术 0.0523。

5.3 最大夏普比率组合

最大夏普比率组合是指单位风险下收益最高的投资组合,是风险与收益权衡的最优解。

# 找出最大夏普比率组合
max_sr <- portfolio_values[which.max(portfolio_values$SharpeRatio), ]

cat("\n=== 最大夏普比率组合 ===\n")
## 
## === 最大夏普比率组合 ===
cat("夏普比率:", round(max_sr$SharpeRatio, 4), "\n")
## 夏普比率: -0.5894
cat("年化收益:", round(max_sr$Return, 4), "\n")
## 年化收益: -0.4156
cat("年化风险:", round(max_sr$Risk, 4), "\n")
## 年化风险: 0.7051
cat("\n权重分配:\n")
## 
## 权重分配:
sr_weights <- round(c(max_sr$Hualing, max_sr$Weichuang, max_sr$Chenguang, max_sr$Lechuang), 4)
names(sr_weights) <- c("华岭股份", "微创光电", "辰光医疗", "乐创技术")
print(sr_weights)
## 华岭股份 微创光电 辰光医疗 乐创技术 
##   0.0280   0.0212   0.0165   0.9343
# 最大夏普比率组合权重直方图
p3 <- data.frame(
  Asset = c("华岭股份", "微创光电", "辰光医疗", "乐创技术"),
  Weights = sr_weights
) %>%
  mutate(Asset = factor(Asset, levels = Asset[order(Weights)])) %>%
  ggplot(aes(x = Asset, y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = '资产', y = '权重', title = "最大夏普比率组合权重") +
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  geom_text(aes(label = paste0(round(Weights * 100, 1), "%")), 
            vjust = -0.5, size = 4) +
  theme(legend.position = "bottom")

ggplotly(p3) %>% layout(title = "最大夏普比率组合权重分配")

最终我们发现最大夏普比率组合的权重分配是:华岭股份 0.0280 微创光电 0.0212 辰光医疗 0.0165 乐创技术 0.9343。

六、有效前沿分析

为了更直观地展示不同投资组合在风险与收益之间的关系,本文基于蒙特卡洛模拟结果绘制有效前沿。

# 创建标签数据框
special_points <- rbind(
  data.frame(Risk = min_var$Risk, Return = min_var$Return, Type = "最小风险组合"),
  data.frame(Risk = max_return$Risk, Return = max_return$Return, Type = "最大收益组合"),
  data.frame(Risk = max_sr$Risk, Return = max_sr$Return, Type = "最大夏普比率")
)

# 使用fill和color分开
p4 <- ggplot() +
  # 所有随机组合 - 使用color映射夏普比率
  geom_point(data = portfolio_values, 
             aes(x = Risk, y = Return, color = SharpeRatio), 
             alpha = 0.5, size = 1.5) +
  scale_color_gradient(low = "blue", high = "red", name = "夏普比率") +
  # 特殊点 - 使用fill映射
  geom_point(data = special_points, 
             aes(x = Risk, y = Return, fill = Type), 
             size = 4, shape = 21, color = "black", stroke = 1.2) +
  scale_fill_manual(values = c("最小风险组合" = "#2ca02c", 
                               "最大收益组合" = "#d62728",
                               "最大夏普比率" = "#ff7f0e"),
                    name = "特殊组合") +
  # 标签
  geom_text(data = special_points, 
            aes(x = Risk, y = Return, label = Type),
            hjust = -0.1, vjust = -0.8, size = 3.5) +
  theme_classic() +
  scale_y_continuous(labels = scales::percent, name = "年化收益率") +
  scale_x_continuous(labels = scales::percent, name = "年化风险") +
  labs(title = "投资组合优化与有效前沿",
       subtitle = paste0("基于", num_port, "个随机组合的蒙特卡洛模拟")) +
  theme(legend.position = "right",
        plot.title = element_text(hjust = 0.5, size = 14),
        plot.subtitle = element_text(hjust = 0.5, size = 10))

# 显示交互式图表
ggplotly(p4) %>% layout(title = "投资组合优化与有效前沿")