本实验中通过四组股票资产的10000次蒙特卡洛模拟，得出固定收益下最低风险与固定风险情形下最大收益，绘制资产有效边界

Downloading data

First lets load our packages

下载数据

加载我们的包

library(PerformanceAnalytics)
## 载入需要的程序包：xts
## 载入需要的程序包：zoo
## 
## 载入程序包：'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## 载入程序包：'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
library(quantmod) 
## 载入需要的程序包：TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(xts) 
library(dygraphs)
library(tibble)
library(tidyr) 
library(ggplot2) 
library(highcharter) 
## Highcharts (www.highcharts.com) is a Highsoft software product which is
## not free for commercial and Governmental use
library(tidyquant)
## 载入需要的程序包：lubridate
## 
## 载入程序包：'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(plotly) # To create interactive charts
## 
## 载入程序包：'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(timetk) # To manipulate the data series
library(forcats)

选择以我学号后四位为列数的资产进行构建，我选择了以下四个资产——

NS(诺思兰德) TH(同辉信息) TR(天润科技) SEAE(万通)

Read the price data.

##选取投资组合对象
#NS(诺思兰德) TH(同辉信息) TR(天润科技) SEAE(万通) 
##读取数据
price <-read.csv("D:/2024_Spring_Teaching/new/price1.csv")
#修改数据格式
price[,1] = as.POSIXct(price[,1])
price = xts(price[,-1], order.by=price[,1])
##绘制走势图
names(price)<-c("NS","TH","TR","SEAE")
dygraph(price)

*接下来，我们将计算这些股票的每日对数回报率。

pr_daily <- to.daily(price, indexAt = "lastof", OHLC = FALSE) 
#计算回报
rt_daily <- na.omit(Return.calculate(pr_daily, method = "log"))
dygraph(rt_daily)

*接下来，让我们计算每种资产的平均每日回报率。计算所有这些股票的协方差矩阵。我们将通过乘以 252 将其年化。

mean_ret <- colMeans(rt_daily)
print(round(mean_ret, 4))
##      NS      TH      TR    SEAE 
## -0.0043 -0.0025 -0.0034 -0.0026
cov_mat <- cov(rt_daily) * 252
print(round(cov_mat,4))
##          NS     TH     TR   SEAE
## NS   0.3079 0.0067 0.0852 0.0575
## TH   0.0067 0.4571 0.0285 0.0151
## TR   0.0852 0.0285 0.3491 0.0773
## SEAE 0.0575 0.0151 0.0773 0.2233

为了计算投资组合的回报和风险（标准差），我们需要

# Calculate the random weights
tick <- c("NS","TH","LECH","SEAE")
wts <- runif(n = length(tick))
wts <- wts/sum(wts)
print(wts)
## [1] 0.2683650 0.2090093 0.2460073 0.2766184
# Calculate the portfolio returns
port_returns <- (sum(wts * mean_ret) + 1)^252 - 1
port_risk <- sqrt(t(wts) %*% (cov_mat %*% wts))
sharpe_ratio <- port_returns/port_risk

我们拥有执行优化所需的一切。我们现在需要做的就是在 20000 个随机投资组合上运行这段代码。为此，我们将使用 for 循环。

在此之前，我们需要创建空向量和矩阵来存储我们的值。

# Creating a matrix to store the weights
num_port <- 20000

# Creating a matrix to store the weights创建矩阵保存权重

all_wts <- matrix(nrow = num_port,
                  ncol = length(tick))

# Creating an empty vector to store创建空向量初始化
# Portfolio returns把它称作投资组合回报

port_returns <- vector('numeric', length = num_port)

# Creating an empty vector to store
# Portfolio Standard deviation

port_risk <- vector('numeric', length = num_port)

# Creating an empty vector to store
# Portfolio Sharpe Ratio

sharpe_ratio <- vector('numeric', length = num_port)

循环%%%%%%%

for (i in seq_along(port_returns)) {
  
  wts <- runif(length(tick))
  wts <- wts/sum(wts)
  
  # Storing weight in the matrix
  all_wts[i,] <- wts
  
  # Portfolio returns
  
  port_ret <- sum(wts * mean_ret)
  port_ret <- ((port_ret + 1)^252) - 1
  
  # Storing Portfolio Returns values
  port_returns[i] <- port_ret
  
  
  # Creating and storing portfolio risk
  port_sd <- sqrt(t(wts) %*% (cov_mat  %*% wts))
  port_risk[i] <- port_sd
  
  # Creating and storing Portfolio Sharpe Ratios
  # Assuming 0% Risk free rate
  
  sr <- port_ret/port_sd
  sharpe_ratio[i] <- sr
}

All the heavy lifting has been done and now we can create a data table to store all the values together. 大功告成！现在创建一个表将所有数据存储在一起

# Storing the values in the table
portfolio_values <- tibble(Return = port_returns,
                           Risk = port_risk,
                           SharpeRatio = sharpe_ratio)


# Converting matrix to a tibble and changing column names
all_wts <- tk_tbl(all_wts)
## Warning in tk_tbl.data.frame(as.data.frame(data), preserve_index, rename_index,
## : Warning: No index to preserve. Object otherwise converted to tibble
## successfully.
colnames(all_wts) <- colnames(pr_daily)
# Combing all the values together
portfolio_values <- tk_tbl(cbind(all_wts, portfolio_values))
## Warning in tk_tbl.data.frame(cbind(all_wts, portfolio_values)): Warning: No
## index to preserve. Object otherwise converted to tibble successfully.

我们有每种资产的权重、风险和回报以及每个投资组合的夏普比率。

接下来，让我们看看最重要的投资组合。

min_var <- portfolio_values[which.min(portfolio_values$Risk),]
max_sr <- portfolio_values[which.max(portfolio_values$SharpeRatio),]
max_rt <- portfolio_values[which.max(portfolio_values$Return),]

绘制每个投资组合的权重。首先是最小方差组合。

p <- min_var %>%
  gather(NS:SEAE, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Assets', y = 'Weights', title = "Minimum Variance Portfolio Weights") +
  scale_y_continuous(labels = scales::percent) 


ggplotly(p)

我试图导出最大回报组合的权重直方图

p <- max_rt %>%
  gather(NS:SEAE, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Assets', y = 'Weights', title = "Maximum Return Portfolio Weights") +
  scale_y_continuous(labels = scales::percent) 


ggplotly(p)

接下来，让我们看一下切线投资组合或夏普比率最高的投资组合。

p <- max_sr %>%
  gather(NS:SEAE, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Assets', y = 'Weights', title = " Tangency Portfolio Weights") +
  scale_y_continuous(labels = scales::percent) 


ggplotly(p)

最后，让我们绘制所有随机投资组合并可视化有效边界。-

p <- portfolio_values %>%
  ggplot(aes(x = Risk, y = Return, color = SharpeRatio)) +
  geom_point() +
  theme_classic() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = 'Annualized Risk',
       y = 'Annualized Returns',
       title = "Portfolio Optimization & Efficient Frontier") +
  #最小方差投资组合点
  geom_point(aes(x = Risk,
                 y = Return), data = min_var, color = 'red') +
  #最大回报率投资组合点
  geom_point(aes(x = Risk,
                 y = Return), data = max_rt, color = 'green') +
  #最大夏普比率投资组合点
  geom_point(aes(x = Risk,
                 y = Return), data = max_sr, color = 'yellow') +
  #添加标签
  annotate('text', x = 0.57, y = -0.44, label = "Tangency Portfolio") +
  annotate('text', x = 0.3, y = -0.49, label = "Minimum variance portfolio") +
  annotate('text', x = 0.43, y = -0.43, label = "Maximum Return portfolio") +
  annotate(geom = 'segment', x = 0.33, xend = 0.33,  y = -0.54 , 
           yend = -0.5, color = 'red', arrow = arrow(type = "open")) +
  annotate(geom = 'segment', x = 0.53, xend = 0.45,  y = -0.47, 
           yend = -0.44, color = 'green', arrow = arrow(type = "open")) +
  annotate(geom = 'segment', x = 0.64, xend = 0.62,  y = -0.47, 
           yend = -0.44, color = 'yellow', arrow = arrow(type = "open")) 
  

ggplotly(p)

后话：当我尝试使用老师的案例R文件生成markdown格式时，我发现我的r语言没有反应，在尝试了检查代码/修改UTF-8格式/搜索文心一言之后（一开始是r版本太低。为此我用了两个小时寻找如何快速的更换包的路径，最后选择了直接用RGUI进行快捷安装，还好没有出现镜像问题），我放弃了抵抗…… 直到现在我也不知道为什么我的文件它不画图 *ps:我把我的微软Office正版化了！