一直想学习一下金融方面的数据挖掘技术，但是对金融理论非常匮乏，所以开始从最基础的“现代金融理论”学起。作为现代金融理论奠基的马科维茨理论当然是必须要好好研究一下的，为了更好的掌握理论，下面就尝试用R进行”均值方差投资组合“有效前沿求解。

用到的知识

• 二次规划（Quadratic Programming）

  二次规划是非线性规划中的一类特殊数学规划问题，在很多方面都有应用，如投资组合、约束最小二乘问题的求解、序列二次规划在非线性优化问题中应用等。在过去的几十年里，二次规划已经成为运筹学、经济数学、管理科学、系统分析和组合优化科学的基本方法。

  二次规划问题可以以下形式来描述：
  min

  受到一个或更多如下型式的限制条件:
  
  

  在R中用quadprog包中的solve.QP来求解二次规划问题。

  函数：solve.QP(Dmat, dvec, Amat, bvec, meq=0, factorized=FALSE)
  详细描述

• 马科维茨投资组合理论
  马科维茨资产组合理论中的最优化资产组合必须符合以下条件之一：

  • 1 在相同预期收益率下，风险最小，即var(ω’r) = ω‘∑ω 达到最小。

  • 2 在相同风险下，预期收益率最大,即max(w’μ)达到最大。
    其中w表示各个风险资产的权重向量
    μ表示各个风险资产的期望收益率向量

    根据上面的条件，最优化资产组合求解可以转化成对min(f(x)) = ω‘∑ω - w’μ 的二次规划求解问题。通过求解min(f(x))可以得到一组风险资产权重向量w。如果给定一个风险因子0<p<1,通过求解 min(f(x)) = ω‘∑ω - pw’μ,我们可以得到多组投资组合，这些投资组合组成的曲线就是投资组合“有效前沿”。

准备数据

使用quantmod包获取股票数据非常方便，直接上代码。
获取5个数据：”腾讯”,”金山”,”中国移动”,”恒生指数”,”招商银行”。

library(quantmod)
stocks = c("0700.hk", "3888.hk","0941.hk","^HSI","3968.hk")
stockNames = c("腾讯","金山","中国移动","恒生指数","招商银行")

getMyStocksData <- function(stocks = c("")){
  stocksdata = list()
  for(i in 1:length(stocks)){
    stock = sprintf('s%s',stocks[i])
    print(paste("beigin get :", stock))
    stocksdata[[stock]] <- getSymbols(stocks[i],from = "2014-06-01",to = Sys.Date(),src = "yahoo",auto.assign=FALSE)
    print(paste("end get :", stock))
  }
  # 对齐日期
  if( length(stocksdata) >= 2){
    ids <- as.character(index(stocksdata[[1]]))
    for( i in 2:length(stocksdata)){
      ids <- intersect(ids, as.character(index(stocksdata[[i]])))
    }
    for(i in 1:length(stocksdata)){
      stock = sprintf('s%s',stocks[i])
      stocksdata[[stock]] <- stocksdata[[stock]][ids]
    }
  }  
  return(stocksdata)
}

getMyStocksReturn <- function(stocksdata){
  for(i in 1:length(stocksdata)){
    ret <- diff(log(stocksdata[[i]][,6]))
    if( i == 1){
      stocksReturn <- data.frame(ret)
    }else{
      stocksReturn <- cbind(stocksReturn, data.frame(ret))
    }
  }
  return(stocksReturn[-1,])
}

sdata <- getMyStocksData(stocks) # stock data
sdata.ret <- getMyStocksReturn(sdata) # stock returns
names(sdata.ret) <- stockNames

执行上面代码，获取到股票收益率数据放到sdata.ret变量里面。查看如下：

> head(sdata.ret)
                   腾讯         金山     中国移动     恒生指数     招商银行
2014-06-04 -0.025317487 -0.002195391 -0.009133588 -0.006000012  0.004101002
2014-06-05 -0.007352879  0.008752791  0.001307992 -0.001817966  0.001363745
2014-06-06 -0.001847162 -0.015368030 -0.009885919 -0.006889210  0.000000000
2014-06-09  0.021938451 -0.006659292  0.004630109  0.007227130  0.004078698
2014-06-10  0.025882656  0.006659292  0.011140788  0.008540040  0.010796124
2014-06-11  0.025230472  0.019715863 -0.021079063 -0.002510088 -0.001342896

求解过程

根据前面的介绍我们知道是要求解min(f(x)) = ω‘∑ω - pw’μ， 0<p<1。通过调整p值，来获取有效前沿组合。

solve.QP 函数原型:

solve.QP(Dmat, dvec, Amat, bvec, meq=0, factorized=FALSE)

设置函数参数：

risk.p <- 0.5
Dmat <- cov(sdata.ret) # 风险资产协方差矩阵
dvec <- colMeans(sdata.ret) * risk.p # 指定风险水平下的期望收益率 
meq <- 1

# 约束条件: sum(x_i) = 1
Amat <- matrix(1, nrow=nrow(Dmat))
bvec <- 1

运行代码:

> ret <- solve.QP(Dmat, dvec, Amat, bvec, meq)
> ret$solution
[1]  2.5103945  0.2184375  3.2731176 -7.0416422  2.0396925

从结果看出，有些权重小于0，而且还有大于1，所以我们还需要再添加一些约束条件：风险资产权重必须大于0且小于1。
修为约束为:

# 约束条件: sum(x_i) = 1 & x_i >= 0
Amat <- cbind(1,diag(nrow(Dmat)))
bvec <- c(1,rep(0,nrow(Dmat)))

运行结果：

> ret <- solve.QP(Dmat, dvec, Amat, bvec, meq)
> ret$solution
[1] 0.0000000 0.0000000 0.2502792 0.0000000 0.7497208

现在这个结果虽然符合上面的约束，但是权重过于集中，有些竟然达到74.9%，对我们来说可能风险过大，所以我们有必要再加一个约束：每个风险资产权重不超过30%比例。

然后我们通过调整p值，获取有效前沿，代码如下：

j <- 0
eq <- vector()
va <- vector()
for( i in seq(0,1.0,0.01)){
  ret <- solve.QP(Dmat, dvec*i, Amat, bvec, meq, factorized=FALSE)
  ret$solution[abs(ret$solution) <= 1e-7] <- 0
  #print(ret$solution)
  #print(paste(t(dvec)%*%ret$solution*252,sqrt(t(ret$solution)%*%Dmat%*%ret$solution*252)))
  eq[j] <- t(dvec)%*%ret$solution*252
  va[j] <- sqrt(t(ret$solution)%*%Dmat%*%ret$solution*252)
  j = j+1
}
eff <- data.frame(eq=eq,va=va)

其中eff就是我们求得的有效前沿在“均值-方差”坐标系中的坐标位置。

绘图展示

利用ggplot绘图，展示eff曲线。

# Color Scheme
ealred  <- "#7D110C"
ealtan  <- "#CDC4B6"
eallighttan <- "#F7F6F0"
ealdark  <- "#423C30"
img <- ggplot(eff, aes(x=va, y=eq)) + geom_point(alpha=.5, color=ealdark) +
  geom_point(data=eff,aes(x=va[eq/va==max(eq/va)],y=eq[eq/va==max(eq/va)]),color=ealred,size=3) +
  geom_line(data=eff,aes(x=va,y=max(eq/va)*va),color=ealtan)+
  ggtitle("Efficient Frontier\nand Optimal Portfolio") + labs(x="Risk (standard deviation of portfolio variance)", y="Return") +
  ggthemes::theme_wsj()
print(img)

图表展示:

完整代码

就到这吧，实在写不下去了。。。

library(quadprog)
library(ggplot2)
library(ggthemes)

############### prepare stocks data begin
library(quantmod)
stocks = c("0700.hk", "3888.hk","0941.hk","^HSI","3968.hk")
stockNames = c("腾讯","金山","中国移动","恒生指数","招商银行")

getMyStocksData <- function(stocks = c("")){
  stocksdata = list()
  for(i in 1:length(stocks)){
    stock = sprintf('s%s',stocks[i])
    print(paste("beigin get :", stock))
    stocksdata[[stock]] <- getSymbols(stocks[i],from = "2014-06-01",to = Sys.Date(),src = "yahoo",auto.assign=FALSE)
    print(paste("end get :", stock))
  }
  # 对齐日期
  if( length(stocksdata) >= 2){
    ids <- as.character(index(stocksdata[[1]]))
    for( i in 2:length(stocksdata)){
      ids <- intersect(ids, as.character(index(stocksdata[[i]])))
    }
    for(i in 1:length(stocksdata)){
      stock = sprintf('s%s',stocks[i])
      stocksdata[[stock]] <- stocksdata[[stock]][ids]
    }
  }  
  return(stocksdata)
}

getMyStocksReturn <- function(stocksdata){
  for(i in 1:length(stocksdata)){
    ret <- diff(log(stocksdata[[i]][,6]))
    if( i == 1){
      stocksReturn <- data.frame(ret)
    }else{
      stocksReturn <- cbind(stocksReturn, data.frame(ret))
    }
  }
  return(stocksReturn[-1,])
}

sdata <- getMyStocksData(stocks) # stock data
sdata.ret <- getMyStocksReturn(sdata) # stock returns
names(sdata.ret) <- stockNames

############### prepare stocks data end


# min(-d^T b + 1/2 b^T D b)
# constraints A^T b >= b_0.
# solve.QP(Dmat, dvec, Amat, bvec, meq=0, factorized=FALSE)

risk.p <- 0.5
Dmat <- cov(sdata.ret)
dvec <- as.vector(colMeans(sdata.ret))
meq <- 1

# Constraints: sum(x_i) = 1
Amat <- matrix(1,nrow=nrow(Dmat))
bvec <- 1

# Constraints: sum(x_i) = 1 & x_i >= 0
Amat <- cbind(1,diag(nrow(Dmat)))
bvec <- c(1,rep(0,nrow(Dmat)))

# Constraints: sum(x_i) = 1 & x_i >= 0 & x_i <= 0.15
#Amat <- cbind(1,diag(nrow(Dmat)),-1 * diag(nrow(Dmat)))
#bvec <- c(1, rep(0, nrow(Dmat)), rep(-0.3, nrow(Dmat)))

j <- 0
eq <- vector()
va <- vector()
for( i in seq(0,1.0,0.01)){
  ret <- solve.QP(Dmat, dvec*i, Amat, bvec, meq, factorized=FALSE)
  ret$solution[abs(ret$solution) <= 1e-7] <- 0
  #print(ret$solution)
  #print(paste(t(dvec)%*%ret$solution*252,sqrt(t(ret$solution)%*%Dmat%*%ret$solution*252)))
  eq[j] <- t(dvec)%*%ret$solution*252
  va[j] <- sqrt(t(ret$solution)%*%Dmat%*%ret$solution*252)
  j = j+1
}
eff <- data.frame(eq=eq,va=va)

# Color Scheme
ealred  <- "#7D110C"
ealtan  <- "#CDC4B6"
eallighttan <- "#F7F6F0"
ealdark  <- "#423C30"
img <- ggplot(eff, aes(x=va, y=eq)) + geom_point(alpha=.5, color=ealdark) +
  geom_point(data=eff,aes(x=va[eq/va==max(eq/va)],y=eq[eq/va==max(eq/va)]),color=ealred,size=3) +
  geom_line(data=eff,aes(x=va,y=max(eq/va)*va),color=ealtan)+
  ggtitle("Efficient Frontier") + 
  labs(x="Risk", y="Return") +
  ggthemes::theme_wsj()
print(img)