本节书摘来异步社区《量化金融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股票的价格，Quandl是一种公共服务，提供了轻松获取大批定量数据的途径。
> 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"):
如果你没有使用认证令牌，可能会出现上面的警告信息。
> 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中。
