本节书摘来异步社区《量化金融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中。