本节书摘来异步社区《量化金融R语言初级教程》一书中的第2章，第2.3节，作者： 【匈牙利】Gergely Daróczi（盖尔盖伊） , 等 译者： 高蓉 , 李茂 责编： 胡俊英，更多章节内容可以访问云栖社区“异步社区”公众号查看。

2.3　使用真实数据

投资组合的最优化完全整合在随后我们会讨论的多个R包中，知道这一点非常有用。但是，跑步之前***先会走路。因此，我们从一个简单的自制R函数开始，我们把它一行行地在下面列出来。

minvariance <- function(assets, mu = 0.005) {
　　 return <- log(tail(assets, -1) / head(assets, -1))
　　 Q　　　<- rbind(cov(return), rep(1, ncol(assets)),
　　　　　　　 colMeans(return))
　　 Q　　　<- cbind(Q, rbind(t(tail(Q, 2)), matrix(0, 2, 2)))
　　 b　　　<- c(rep(0, ncol(assets)), 1, mu)
　　 solve(Q, b)
}

这就是我们在“（拉格朗日）定理”这一节讨论的算法的直接运用。

为了做出示范，我们从Quandl超集中选取了一些IT股票的价格(http://www.quandl.com/USER_1KR/1KT)，Quandl是一种公共服务，提供了轻松获取大批定量数据的途径。这个超链接［（http://www.quandl.com/api/v1/datasets/USER_1KR/1KT.csv）］指向一个可下载的逗号分隔值（CSV）文件，可以保存在硬盘上，并通过read.csv导入到R中，此外还有一种更直观的方法来完成这种操作，就是借助于包括在上述超链接中的键的帮助。

> library(Quandl)
> IT <- Quandl('USER_1KR/1KT',
+　　　　　 start_date = '2008-01-01', end_date = '2012-12-31')
Warning message:
In Quandl("USER_1KR/1KT", start_date = "2008-01-01", end_date = 
"2012-12-31"):

如果你没有使用认证令牌，可能会出现上面的警告信息。请访问http://www.quandl.com/help/r，或者每天仅从Quandl下载10个数据集。

> str(IT)
'data.frame': 1259 obs. of 6 variables:
$ Date: Date, format: "2008-01-02" "2008-01-03" ...
$ AAPL: num 199 195 191 181 180 ...
$ GOOG: num 693 685 680 654 653 ...
$ MSFT: num 35.8 35.2 35.2 34.5 34.7 ...
$ IBM : num 109 105 104 100 100 ...
$ T　 : num 41.5 41.2 41 41.1 41.3 ...

因此，我们载入Quandl包，它提供的Quandl函数需要以下这几个参数。

第一个参数（code=USER_1KR/1KT）是数据集在Quandl上的代码。
start_date和end_date参数可以任意指定我们感兴趣的时间段，在这里我们设置为从五年前到现在。
更多的选择请参考?Quandl。例如，可以使用type导入已经存储在某些时间序列对象中的数据，来代替最初的data.frame。
在新创建的IT变量上运行str命令，结果显示了R对象的内部结构，它当前包括一个Date域和5个数值格式的资产价格。

将价格从IT中指派给assets之后（不包括第一列Date），我们来逐行运行之前的minvariance函数的主体部分。首先，我们计算资产收益率，除了第一个值，每个值都除以自身的前一个值，再对每个商计算log。

> assets <- IT[, -1]
> return <- log(tail(assets, -1) / head(assets, -1))

请注意，收益率也可以直接由timeSeries包的returns函数计算，为了教学的目标，我们在这里不调用它。为了确认命令做了什么，我们来检查新创建的变量的前一些值。

> head(return)
　　　　 AAPL　　　　GOOG　　　　　MSFT　　　　　IBM　　　　　　T
2 -0.019560774 -0.011044063 -0.0160544217 -0.038916144 -0.0072534167
3 -0.020473237 -0.008161516 -0.0008521517 -0.008429976 -0.0043774389
4 -0.054749384 -0.038621208 -0.0183544011 -0.036242948  0.0007309051
5 -0.006142967 -0.001438475  0.0046202797 -0.001997005  0.0051014322
6 -0.050317921 -0.035793820 -0.0396702524 -0.023154566 -0.0514590970
7  0.036004806  0.023482511  0.0292444412 -0.003791959 -0.0123204844

接着，我们来建立拉格朗日定理指定的线性等式系统的左侧：left[ {begin{array}{*{20}{c}}Q{overrightarrow 1 }\{overrightarrow r }end{array}} right]，其中我们按行（rbind）把协方差矩阵（cov），根据数据集列数（ncol）重复（rep）次数的1的组合向量，以及收益率的均值（colMeans）结合在一起。

> Q <- rbind(cov(return), rep(1, ncol(assets)), colMeans(return))

最后，结果如下。

> round(Q, 5)
　　　　AAPL　　GOOG　　MSFT　　 IBM　　　　 T
AAPL 0.00063 0.00034　0.00025 0.00023　0.00022
GOOG 0.00034 0.00046　0.00023 0.00019　0.00018
MSFT 0.00025 0.00023　0.00034 0.00018　0.00018
IBM　0.00023 0.00019　0.00018 0.00024　0.00016
T　　0.00022 0.00018　0.00018 0.00016　0.00028
　　 1.00000 1.00000　1.00000 1.00000　1.00000
　　 0.00075 0.00001 -0.00024 0.00044 -0.00018

请注意，为了结果易读，我们已经将运算结果保留至5位小数。还有请注意，微软（Microsoft，MSFT）和美国电话电报公司（AT&T）的均值都为负。现在，我们也需要将矩阵tail的最后两行结合在一起生成新列放在矩阵左部，为了拉格朗日定理（2×2矩阵）指定的线性系统的完整性，再将其他部分设置为零。

> Q <- cbind(Q, rbind(t(tail(Q, 2)), matrix(0, 2, 2)))
> round(Q, 5)
　　　　AAPL　　GOOG　　 MSFT　　 IBM　　　　T
AAPL 0.00063 0.00034　0.00025 0.00023　0.00022 1　0.00075
GOOG 0.00034 0.00046　0.00023 0.00019　0.00018 1　0.00001
MSFT 0.00025 0.00023　0.00034 0.00018　0.00018 1 -0.00024
IBM　0.00023 0.00019　0.00018 0.00024　0.00016 1　0.00044
T　　0.00022 0.00018　0.00018 0.00016　0.00028 1 -0.00018
　　 1.00000 1.00000　1.00000 1.00000　1.00000 0　0.00000
　　 0.00075 0.00001 -0.00024 0.00044 -0.00018 0　0.00000

在默认情形下，mu是0.005（在最小方差函数的参数中指定）。这是线性系统右侧向量left[ {begin{array}{*{20}{c}}0\1\mu end{array}} right]的最后一个值。

> mu <- 0.005
> b <- c(rep(0, ncol(assets)), 1, mu)
> b
[1] 0.000 0.000 0.000 0.000 0.000 1.000 0.005

在成功建立线性系统的各个部分之后，你只剩下了求解的任务。

> solve(Q, b)
　　　　AAPL　　 GOOG　　　　　　　MSFT　　　　　IBM　　　　　　T
2.3130600636  -1.0928257246 -2.7830264740 4.9871895547 -2.4243974197 
-0.0001254637 -1.2082468413

上面的代码等价于一步运行这个函数，它会取数据集和最优的可选收益率作为参数。为了在想要的预期收益率水平下获得最小方差，计算结果给出了最优权重和朗格朗日乘子。

> minvariance(IT[, -1])
　　　　AAPL　　　　　GOOG　　　　 MSFT　　　　 IBM　　　　　　 T
2.3130600636　-1.0928257246 -2.7830264740 4.9871895547 -2.4243974197
-0.0001254637 -1.2082468413

注意，在Microsoft和AT&T的股票前面，Google的最优化头寸为做空。我们可以使用这个结果进而寻找最优化问题的完整解，借助其他软件和“write.csv”函数，这个求解过程还可以进一步深化。除了在给定的收益率水平上计算最小方差，我们还可以对更大范围的收益率求解最小方差。如果我们修正这份代码，可以得到如下的结果。

frontier <- function(assets) {
　　return <- log(tail(assets, -1) / head(assets, -1))
　　Q <- cov(return)
　　n <- ncol(assets)
　　r <- colMeans(return)
　　Q1 <- rbind(Q, rep(1, n), r)
　　Q1 <- cbind(Q1, rbind(t(tail(Q1, 2)), matrix(0, 2, 2)))
　　rbase <- seq(min(r), max(r), length = 100)
　　s <- sapply(rbase, function(x) {
　　　　y <- head(solve(Q1, c(rep(0, n), 1, x)), n)
　　　　y %*% Q %*% y
　　})
　　plot(s, rbase, xlab = 'Return', ylab = 'Variance')
}

除了在资产的最小收益率和最大收益率之间（seq）取过一个不同的收益率值（length = 100），和除了计算了最优投资组合的方差，这份代码与之前的版本一样。因此，我们可以画出收益率-方差对（s和rbase）来表示问题的解。结果显示在图2-2中。

在方差-收益率平面上，理想的收益率-最小方差曲线叫作投资组合前沿（Portfolio Frontier）。忽略它向下方倾斜的部分（同样的方差可以用更高的预期收益率达到），我们得到了有效前沿（Efficient Frontier），毫无疑问必须选择有效前沿上的组合。

众所周知，两个给定的收益率水平就足以计算投资组合前沿，把得到的投资组合连接起来就得到了整个前沿。

相似的结果可以通过R包中的一些内置函数来计算，无需太多代码。例如，fPortfolio包提供了一组有用方法，适合于timeSeries对象的应用问题。为了完成这个目的，我们必须把最初的IT数据集资产列转换为一个根据第一列定义的timeSeries对象。

> library(timeSeries)
> IT <- timeSeries(IT[, 2:6], IT[, 1])

就像我们在均值-方差函数中所做的那样，通过对每一个元素除以前一个值再计算对数，收益率可以定义为时间序列，但是，一些有用的时间序列命令（比如lag）可以更轻松地实现这个目的。

> log(lag(IT) / IT)

或者，使用其他内置函数可以是更简单的方式。

> IT_return <- returns(IT)

现在，我们已经有了一个时间序列对象，因此绘出收益率很容易。

> chart.CumReturns(IT_return, legend.loc = 'topleft', main ='')

IT_return中5个股票的收益率看上去和图2-3所示的一样。

通过绘出portfolioFrontier的结果，可以通过交互方式画出上述的前沿图像（图2-3）。

> library(fPortfolio)
> plot(portfolioFrontier(IT_return))

Make a plot selection (or 0 to exit):

1:　 Plot Efficient Frontier
2:　 Add Minimum Risk Portfolio
3:　 Add Tangency Portfolio
4:　 Add Risk/Return of Single Assets
5:　 Add Equal Weights Portfolio
6:　 Add Two Asset Frontiers [LongOnly Only]
7:　 Add Monte Carlo Portfolios
8:　 Add Sharpe Ratio [Markowitz PF Only]

为了模仿我们在上述代码中实现的内容，我们绘出带有卖空约束的前沿图像。

> Spec = portfolioSpec()
> setSolver(Spec) = "solveRshortExact"
> Frontier <- portfolioFrontier(as.timeSeries(IT_return),
+　　　　　　　　 Spec, > constraints = "Short")
> frontierPlot(Frontier, col = rep('orange', 2), pch = 19)
> monteCarloPoints(Frontier, mcSteps = 1000, cex = 0.25, pch =19)
> grid()

在上面的这段代码中，我们通过一个使无限制卖空卖出的投资组合最优化的函数（solveRshortExact），设置了一个特定的portfolioSpecS4对象。通过橙色圆点（pch = 19）表示的前沿图形（frontierPlot）给出了计算的结果（portfolioFrontier）。一些小（cex =0.25）的蒙特卡罗模拟的点以及作为背景的网格线也添加到了图上，这些都显示在图2-4中。