R求解投资组合有效前沿

一直想学习一下金融方面的数据挖掘技术,但是对金融理论非常匮乏,所以开始从最基础的“现代金融理论”学起。作为现代金融理论奠基的马科维茨理论当然是必须要好好研究一下的,为了更好的掌握理论,下面就尝试用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)

(转)R书精选16本

收藏一下,后面慢慢学习。

文章转自:http://xccds1977.blogspot.com/2013/02/r.html

以前人的烦恼是没有书可读,现在人的烦恼是书太多了。关于R语言的书已经出版很多了,博主大约读过其中的四十多本,但是书在精,而不在多,学在透,而不在速。把有限的时间放到无限的书海中,这不是阅读的真意。本着造福学习者的角度,博主精选出十二本R书。什么是好书的标准?我以为是:有案例,有代码,有习题,有讲解,逻辑清楚,排版精良,体系完备,互有补充,内容千锤百炼,值得反复揣摩。书单均为英文版,都可以从网上找到。当然这份书单的选择是有主观偏见的。

一、初学入门:
《R in Action》
《The Art of_R Programming》
入门者可首选两本,前者从统计角度入手,分高中低三部分由浅入深的讲解了如何用R来实现统计分析,另外此书已经有中文版面世。后者从程序编写的角度入手,对R的本身特点进行了清晰的介绍。中文版应该快有了。
更新:《learning R》
这本书没有单纯的讲语法,而是和数据分析的流程结合了起来,从数据获取到数据整理再到分析和报告,有一气呵成的感觉,此外最后两章讲如何写稳健的R代码以及写包都是非常精彩的。

二、统计进阶:
《A Handbook of Statistical Analyses_Using_R》
《Modern Applied Statistics With S》
这两本书基本上涵盖了统计的一些高阶内容,例如多元分析、多层回归模型、荟萃分析、生存分析等内容。案例丰富,公式不多,值得反复学习参考。

三、科学计算:
《Introduction to Scientific Programming and Simulation Using R》
除了统计分析外,此书独特之处在于使用R来做数值分析,如求根,最优化,数值积分。还包括了一些常见的模拟技术。书后的习题和最后的案例非常有用。该书的中文版据说还在翻译。

四、数据挖掘:
《Data Mining with R_ Learning with Case Studies》
《Machine Learning for Hackers》*
两本侧重于数据挖掘的R书,全是以案例为线索,示范的代码量很大。跟一遍下来会有很大的收获。
更新:
《An Introduction to Statistical Learning》**
这本书可以说是另一本数据挖掘大作《The Elements of Statistical Learning》的R实现手册,体系结构基本一致,更强调用R来实现,更难得的地方是提供了很好的习题。

五、数据绘图:
《ggplot2 Elegant Graphics for Data Analysis》
ggplot2还有什么好说的呢,R中最优秀的绘图包,但由于近期该包升级很快,这书显得有些过时。好在中文版进行了大幅更新,即将面世。
更新:《R Graphics Cookbook》
这本书也是RStudio公司的人出的,似乎是Hadley的学生吧,主要是各种ggplot2包的例子,也包括了用其它包来画图,建议通读一遍。

六、参考手册:
《R Cookbook》
《R in a Nutshell》
有时候我们需要类似词典的案头参考手册,以方便随时查阅。又或者可以通读一遍以查漏补缺。上面两本书虽然有些厚度,但仍然推荐之。后者的中文版也在翻译状态。

七、高级编程:
《R Programming for Bioinformatics》
《software for data analysis programming with R》
如果你是初学者,不要去看上面两本书。如果你想进阶为专家级R用户,那你需要精读它们。前者讲解了R少为人知的一面,例如字符处理、正则表达和XML,还有报错处理以及与其它语言的交互。后者更是编写生产级代码的圣经指南。
更新:《Advanced R programming》
Hadley的力作,只是还没有写完,已经可以从这里看到部分,清楚的讲解了R的函数式编程思想和写R包的各种细节,要迈入R高手,不得不读。

分类器哪家强?

看看别人的评测报告。

Do we Need Hundreds of Classifiers to Solve Real World Classification Problems?

Manuel Fernández-Delgado, Eva Cernadas, Senén Barro, Dinani Amorim; 15(Oct):3133−3181, 2014.

Abstract

We evaluate 179 classifiers arising from 17 families (discriminant analysis, Bayesian, neural networks, support vector machines, decision trees, rule-based classifiers, boosting, bagging, stacking, random forests and other ensembles, generalized linear models, nearest-neighbors, partial least squares and principal component regression, logistic and multinomial regression, multiple adaptive regression splines and other methods), implemented in Weka, R (with and without the caret package), C and Matlab, including all the relevant classifiers available today. We use 121 data sets, which represent the whole UCI data base (excluding the large- scale problems) and other own real problems, in order to achieve significant conclusions about the classifier behavior, not dependent on the data set collection.

The classifiers most likely to be the bests are the random forest (RF) versions , the best of which (implemented in R and accessed via caret) achieves 94.1% of the maximum accuracy overcoming 90% in the 84.3% of the data sets. However, the difference is not statistically significant with the second best, the SVM with Gaussian kernel implemented in C using LibSVM, which achieves 92.3% of the maximum accuracy. A few models are clearly better than the remaining ones: random forest, SVM with Gaussian and polynomial kernels, extreme learning machine with Gaussian kernel, C5.0 and avNNet (a committee of multi-layer perceptrons implemented in R with the caret package). The random forest is clearly the best family of classifiers (3 out of 5 bests classifiers are RF), followed by SVM (4 classifiers in the top-10), neural networks and boosting ensembles (5 and 3 members in the top-20, respectively).

查看报告pdf: http://jmlr.org/papers/volume15/delgado14a/delgado14a.pdf

Python 学习资源

最近写一个基于python的web应用,趁机补补python知识,记录一下收集到的python学习资源(不定期更新)。