Downloading data

First lets load our packages

library(xts)#处理时间序列
## Warning: 程辑包'xts'是用R版本4.2.2 来建造的
## 载入需要的程辑包：zoo
## Warning: 程辑包'zoo'是用R版本4.2.3 来建造的
## 
## 载入程辑包：'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(dygraphs)#绘图
## Warning: 程辑包'dygraphs'是用R版本4.2.3 来建造的
library(tibble)#数据形式
library(ggplot2)#绘图
## Warning: 程辑包'ggplot2'是用R版本4.2.3 来建造的
library(PerformanceAnalytics)#计算收益率与协方差矩阵
## Warning: 程辑包'PerformanceAnalytics'是用R版本4.2.3 来建造的
## 
## 载入程辑包：'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend

Next lets select a few stocks to build our portfolios. We will choose the following 4 stocks TH(诺斯兰德) Wc(齐鲁华信) Sz(同辉信息) HN(万通液压) Read the price data.

#读入数据
price <-read.csv("C:/Users/86133/Desktop/portfolio/data.csv")
#处理时间
t<-as.POSIXct(price[4:313,1])
#把价格由文本转化成数，并且保存成dataframe,将dataframe在转成时间序列
a1<-as.numeric(price[4:313,2])
a2<-as.numeric(price[4:313,19])
a4<-as.numeric(price[4:313,3])
a9<-as.numeric(price[4:313,20])
names<-c("TH","Wc","Sz","HN")
mydata1<-data.frame('TH'=a1,'Wc'=a2,'Sz'=a4,'HN'=a9)
price2= xts(mydata1, order.by=t)
#绘制价格走势
dygraph(price2)

Next we will calculate the daily returns for these stocks. We will use the logarithmic returns.

head(price2)
##               TH   Wc    Sz    HN
## 2022-01-04 25.47 8.26 10.03  9.48
## 2022-01-05 24.16 8.35  9.60 10.30
## 2022-01-06 24.53 8.40  8.50  9.32
## 2022-01-07 23.19 8.32  8.29  9.00
## 2022-01-10 22.30 8.40  8.76  9.03
## 2022-01-11 20.50 8.36  8.30  8.91
pr_daily <- to.daily(price2, indexAt = "lastof", OHLC = FALSE) 
head(pr_daily)
##               TH   Wc    Sz    HN
## 2022-01-31 25.47 8.26 10.03  9.48
## 2022-01-31 24.16 8.35  9.60 10.30
## 2022-01-31 24.53 8.40  8.50  9.32
## 2022-01-31 23.19 8.32  8.29  9.00
## 2022-01-31 22.30 8.40  8.76  9.03
## 2022-01-31 20.50 8.36  8.30  8.91
rt_daily <- na.omit(Return.calculate(pr_daily, method = "log")) 
head(rt_daily)
##                     TH           Wc          Sz           HN
## 2022-01-31 -0.05280292  0.010836951 -0.04381750  0.082959579
## 2022-01-31  0.01519849  0.005970167 -0.12169693 -0.099981267
## 2022-01-31 -0.05617571 -0.009569451 -0.02501619 -0.034938051
## 2022-01-31 -0.03913447  0.009569451  0.05514594  0.003327790
## 2022-01-31 -0.08416179 -0.004773279 -0.05394039 -0.013378126
## 2022-01-31  0.02600050 -0.005998818  0.04706751 -0.004499445
dygraph(rt_daily)

Next let’s calculate the mean daily returns for each asset. Calculate the covariance matrix for all these stocks. We will annualize it by multiplying by 252.

mean_ret <- colMeans(rt_daily)
print(round(mean_ret, 4))
##      TH      Wc      Sz      HN 
## -0.0023 -0.0013 -0.0043 -0.0013
cov_mat <- cov(rt_daily) * 252
print(round(cov_mat,4))
##        TH     Wc     Sz     HN
## TH 0.2788 0.0447 0.1010 0.0432
## Wc 0.0447 0.0970 0.0480 0.0389
## Sz 0.1010 0.0480 0.3079 0.0638
## HN 0.0432 0.0389 0.0638 0.0824

Simulation

er_p<-c()
sigma_p<-c()
sharpe_p<-c()
all_wts <- matrix(nrow = 1000,
                  ncol = 4)
for(i in 1:1000){w<-runif(4,min=0,max=1)
                 w<-w/sum(w)
                 all_wts[i,] <- w
                 miu<-((t(w)%*%mean_ret+1)^252)-1
                 std<-sqrt(t(w)%*%cov_mat%*%w)
                 sharpe_ration<-miu/std
                 sharpe_p<-append(sharpe_p,sharpe_ration,length(sharpe_p))
                 er_p<-append(er_p,miu,length(er_p))
                 sigma_p<-append(sigma_p,std,length(sigma_p))}

find the maxmium-return portfolio,min-variation portfolio and maxmium-Sharperation portfolio

#最值查找
library(tibble)
portfolio_values <- tibble(Return = er_p,
                           Risk = sigma_p,
                           SharpeRatio = sharpe_p,weights=all_wts)
min_var <- portfolio_values[which.min(portfolio_values$Risk),]
max_return<-portfolio_values[which.max(portfolio_values$Return),]
max_sr <- portfolio_values[which.max(portfolio_values$SharpeRatio),]

draw the histogram of the four portfolios we found just now

#最值权重分配直方图
p_return<-data.frame('name'=names,'weights'=as.vector(as.matrix(max_return[,4])))
p_1<-ggplot(data=p_return,aes(x=name,y=weights))+
  geom_bar(stat='identity',fill='pink')+
  theme_minimal() +
  labs( title = "Max ruturn Portfolio Weights")+
  scale_y_continuous(labels = scales::percent) 
print(p_1)

p_var<-data.frame('name'=names,'weights'=as.vector(as.matrix(min_var[,4])))
p_2<-ggplot(data=p_var,aes(x=name,y=weights))+
  geom_bar(stat='identity',fill='pink')+
  theme_minimal() +
  labs( title = "min risk Portifolio Weights")+
  scale_y_continuous(labels = scales::percent) 
print(p_2)

p_sr<-data.frame('name'=names,'weights'=as.vector(as.matrix(max_sr[,4])))
p_3<-ggplot(data=p_sr,aes(x=name,y=weights))+
  geom_bar(stat='identity',fill='pink')+
  theme_minimal() +
  labs( title = "Tangency Portfolio Weights")+
  scale_y_continuous(labels = scales::percent) 
print(p_3)

draw the frontier and mark the four potifolios

#最小方差前沿
mydata2<-data.frame('expected'=er_p,'volatility'=sigma_p,'sharperation'=sharpe_p)
p<-ggplot(data=mydata2,aes(x=volatility,y=expected,color=sharperation))+
  geom_point()+
  geom_point(data=min_var,aes(x=Risk,y=Return),col="red")+
  geom_point(data=max_return,aes(x=Risk,y=Return),col="red")+
  geom_point(data=max_sr,aes(x=Risk,y=Return),col="red")+
  labs(title="最小方差前沿",caption="投资组合")+
  theme_classic()
print(p)