量化金融R语言初级教程

版权信息

书名：量化金融R语言初级教程

ISBN：978-7-115-45123-1

本书由人民邮电出版社发行数字版。版权所有，侵权必究。

您购买的人民邮电出版社电子书仅供您个人使用，未经授权，不得以任何方式复制和传播本书内容。

我们愿意相信读者具有这样的良知和觉悟，与我们共同保护知识产权。

如果购买者有侵权行为，我们可能对该用户实施包括但不限于关闭该帐号等维权措施，并可能追究法律责任。

•　著　　　　[匈牙利] Gergely Daróczi 等

　 译　　　　高　蓉　 李　茂

　 责任编辑　胡俊英

•　人民邮电出版社出版发行　　北京市丰台区成寿寺路11号

　 邮编　100164 　电子邮件　315@ptpress.com.cn

　 网址　http://www.ptpress.com.cn

•　读者服务热线：(010)81055410

　 反盗版热线：(010)81055315

版权声明

Copyright ©2013 Packt Publishing. First published in the English language under the title Introduction to R for Quantitative Finance.

All rights reserved.

本书由英国Packt Publishing公司授权人民邮电出版社出版。未经出版者书面许可，对本书的任何部分不得以任何方式或任何手段复制和传播。

版权所有，侵权必究。

内容提要

R是用于统计分析、绘图的语言和操作环境。它是属于GNU系统的一个自由、免费、源代码开放的软件，是一个用于统计计算和统计制图的优秀工具。

本书通过9章的内容向读者详细介绍使用R语言实现量化金融的一些基础知识和方法，内容包括时间序列分析、投资组合优化、资产定价模型、固定收益证券、估计利率期限结构、衍生品定价、信用风险管理、极值理论和金融网络等。

本书的目标读者是那些希望通过R语言来解决量化金融问题的读者，如果读者具备一定的金融知识，将会对本书的阅读有较大的帮助。通过阅读本书，读者将学习到有关R语言的诸多核心内容，并了解R语言在量化金融方面的各类应用。

译者序

自20世纪50年代马可维茨开创投资组合理论以来，金融学开始与数学模型和数据分析紧密融合，逐渐形成了高度数理化的现代金融学。时至今日，无论是对传统的金融学理论，还是对如日中天的量化投资，又或是对方兴未艾的金融科技（Fintech），数学模型与数据科学源源不断地为这些金融理论和实践提供着发展工具，一起构建了复杂深奥、丰富多彩的量化金融世界。由于数学模型与数据分析本身的复杂性，量化金融技术一直被公众视为数学家或华尔街的“火箭科学家”的秘密武器。

近年来，开源软件得到了迅速推广，促进了各行各业对数据科学的学习与应用。R就是这样一种开源数据分析工具。它灵活易用，数据分析能力强大，从技术上降低了量化金融技术的学习和应用难度，目前在世界范围内被广泛地采用，在中国也拥有大量用户。但是，可以结合量化金融与数据科学的图书，在市场上依然很少，而本套书正是基于这种需求而作的关于量化金融的学习教程。本套书选择R语言作为工具，通过R语言学习量化金融，帮助读者快速掌握知识并加以实践应用，为各种金融问题提供实践解决方案，同时还可使用第三方贡献的免费R包。

本套书内容广泛、取舍精当，涵盖了实证金融、金融工程、交易策略与银行管理等内容，可以帮助读者全面学习金融理论与技术。本套书既包含金融时间序列分析、资产定价和期权定价等理论知识，也包括投资组合管理、信用风险管理等实践知识，还涉及金融学的前沿领域以及2008年金融危机后发展起来的金融网络分析理论。本套书包括两本，分别是《量化金融R语言初级教程》和《量化金融R语言高级教程》。学习《量化金融R语言初级教程》不需要读者精通R和金融理论，但需要对金融领域具有一定了解。学习《量化金融R语言高级教程》要求读者具备一定的R编程能力和金融学基本概念，建议在学习《量化金融R语言初级教程》的基础上再学习《量化金融R语言高级教程》。

作为金融学教师，我们认为本套书是通向量化金融世界的一把“金钥匙”，可以有效促进金融知识和计算技术的传播。我们很荣幸地承担了本书的翻译工作。在此特别感谢胡俊英编辑，胡编辑认真专业的审核工作，有力保证了本书翻译的如期完成，我们亦受益匪浅。同时特别感谢天津理工大学的李茂老师，书中多处的专业词汇敲定，都来自与李茂老师的讨论结果。当然，文责自负。同时，感谢杭州电子科技大学2016年高等教育研究资助项目YB201631“投资学教学与R软件应用”的支持。

——高蓉，2017年3月

译者简介

高蓉，任教于杭州电子科技大学经济学院。博士毕业于南开大学经济学院，本科毕业于南开大学数学科学学院。研究领域包括资产定价、实证金融、数据科学应用。已出版教材《实验投资学》，译著《数据科学入门》，发表学术论文数篇。

李茂，任教于天津理工大学数学系，毕业于北京师范大学。热爱数据科学，从事与统计和数据分析相关的教学和研究工作。

作者简介

Gergely Daróczi是一位社会学博士学位候选人，拥有大约8年R编程的数据管理和分析任务的工作经验。Gergely数年来在匈牙利的多所大学讲授统计学课程并从事数据分析工作，最近他还创建并协调着一个总部位于英国的在线报告创业公司。后者作为一种服务性平台，其软件或者平台称为rapporter.net，可以对本书涉及的所有方法和技术提供一个直观的界面和接口。他对本书的贡献是提供了量化金融问题和方法的R实现。

“我非常感谢家庭成员的支持和理解。尽管为了完成本书的工作，他们已经很久没有见到我，但他们也依然默默支持着我。我还衷心感谢所有在匈牙利布达佩斯考文纽斯大学任教的诸位合作者，他们为本次合作提供了有用的内容。”

Michael Puhle在德国帕绍大学（University of Passau）获得了金融学博士学位。他曾在慕尼黑的安联资产管理公司（Allianz Global Investors）担任高级风险控制经理多年，后来在毕马威金融风险管理（KPMG's Financial Risk Management）部门中担任管理助理，在那里他就市场风险模型为银行提供咨询。他还是斯普林格出版社（Springer Publishing）出版的《债券组合优化》（Bond Portfolio Optimization）的作者之一。

Edina Berlinger是毕业于布达佩斯考文纽斯大学（Corvinus University of Budapest）的经济学博士。她是一名助理教授，讲授公司金融、投资学和金融风险管理。她还担任大学金融系的领导职务，也是匈牙利科学院金融分委员会的主席。她的专业涉及学生贷款系统、风险管理，最近又涉及了网络分析领域。在学生贷款设计、流动性管理、异质代理模型和系统风险方面，她领导过一些研究项目。

Péter Csóka是布达佩斯考文纽斯大学的助理教授，同时也是匈牙利科学院经济与区域研究中心、博弈论研究组的研究员。2008年，他在马斯特里赫特大学（Maastricht University）获得了博士学位。他的研究主题包括风险管理、风险资本配置、博弈论、公司金融以及一般均衡理论。他目前致力于分析对系统风险和非流动性资产组合的风险贡献。他在《运筹研究的数学方法》（Mathematical Methods of Operational Research）、《欧洲运筹研究杂志》（European Journal of Operational Research）、《博弈与经济行为》（Games and Economic Behaviour）以及《银行与金融杂志》（Journal of Banking and Finance）发表过论文。他还是布达佩斯金融市场流动性年度会议组织委员会的主席。

Daniel Havran是匈牙利科学院经济研究所、经济与区域研究中心的博士后研究员。他还在布达佩斯考文纽斯大学担任兼职的助理教授，在那里，他讲授公司金融（本科和博士水平）和信用风险管理（硕士水平）的课程。2011年，他在布达佩斯考文纽斯大学获得了经济学博士学位。他研究的兴趣方向是公司现金、基金流动性管理以及场外市场的信用衍生品。

Márton Michaletzky在2011年从布达佩斯考文纽斯大学获得了经济学博士学位。在2000～2003年间，他是协和证券有限公司（Concorde Securities Ltd）的风险经理和宏观经济分析师。作为资本市场交易经理，他在匈牙利国家高速公路管理公司获得了30亿欧元的证券化经验。2012年，他参与了一次IPO的准备工作以及匈牙利金融服务提供商的私人配售。在加入DBH投资之前，他是CUB金融系的助理教授。

Zsolt Tulassay在一家专业的美国投资银行任量化分析师，从事评估衍生品定价模型相关的工作。在此之前，Zsolt是布达佩斯考文纽斯大学金融系的助理讲师，讲授衍生品、量化风险管理和金融计量经济学。Zsolt拥有布达佩斯考文纽斯大学和中欧大学的硕士学位。他研究的兴趣方向包括衍生品定价、构建收益率曲线、流动性风险以及异质代理模型。

Kata Váradi自2013年以来一直在布达佩斯考文纽斯大学任金融学的助理教授。2009年，Kata从布达佩斯考文纽斯大学研究生毕业，并于2012年获得了博士学位，她论文的主题是关于匈牙利股票市场的市场流动性风险分析。她的研究领域包括市场流动性、固定收益证券以及医疗系统的网络。除了研究，她对教学也很积极。主要讲授公司金融、投资学、估值以及跨国金融管理。

Agnes Vidovics-Dancs是博士学位候选人和布达佩斯考文纽斯大学金融系的助理教授。在此之前，她是匈牙利政府债务管理局的初级风险经理。她的主要研究领域是通常的政府债务管理，特别是主权危机和违约。

审稿人简介

Hari Shanker Gupta博士是一位算法交易系统开发领域的高级量化研究分析师。在此之前，他是位于印度班加罗尔的印度科学研究院（Indian Institute of Science，IISc）的博士后。在2010年，他在IISc获得了应用数学和科学计算领域的博士学位，他还在印度瓦拉纳西的巴纳拉斯印度大学（Banaras Hindu University，BHU）完成了硕士学位。在瓦拉纳西，因为其杰出的表现，他获得过4次金奖。

Hari已经在数学和科学计算的著名期刊上发表过5篇研究论文。他拥有在数学、统计和计算领域的工作经验，涉及数值方法、偏微分方程、数理金融、随机分析、数据分析、时间序列分析、有限差分以及有限元方法等领域。他擅长数学软件Matlab、统计编程语言R、编程语言C，目前，他在Python平台下工作。

Ronald Hochreiter是维也纳大学经贸大学的金融、会计和统计系的助理教授。在2005年，他在维也纳大学获得了他的计算管理科学的博士学位。他是一个热心的R用户，并开发了主要用于最优化建模和金融应用的R包。他的R项目汇总可以在http://www.hochreiter.net/R/找到，而且他的一些关于R使用的金融工程在线教程可以在http://www.finance-r.com/找到。

前言

本书将向你讲述如何使用统计计算语言R和量化金融知识来解决真实世界的量化金融问题。本书包括了丰富的主题，从时间序列分析到金融网络。每章都会简要地介绍理论知识并使用R来解决一个具体问题。

本书内容

第1章“时间序列分析”（Michael Puhle），介绍了用R处理时间序列数据。并且，你会学到如何建模和预测房价，使用协整改善对冲比，以及对波动率建模。

第2章“投资组合优化”（Péter Csóka，Ferenc Illés，Gergely Daróczi），包括了投资组合选择背后的理论思想，并说明了如何将这些知识运用于真实世界的数据。

第3章“资产定价模型”（Kata Váradi，Barbara Mária Dömötör，Gergely Daróczi），建立在前一章的基础上，给出了刻画资产收益率与风险之间关系的模型。本章包括资本资产定价模型和套利定价理论。

第4章“固定收益证券”（Márton Michaletzky，Gergely Daróczi），是处理固定收益产品的基础。在这一章中，你会学到如何计算这类产品的风险，以及构建对利率变化免疫的投资组合。

第5章“估计利率期限结构”（Tamás Makara，Gergely Daróczi），介绍了收益率曲线的概念，并说明了如何使用政府债券的价格估计收益率曲线。

第6章“衍生品定价”（Ágnes Vidovics-Dancs，Gergely Daróczi），使用离散和连续时间模型解释了衍生品定价。而且，你还会学到如何计算衍生品风险的度量以及所谓的“希腊字母”。

第7章“信用风险管理”（Dániel Havran，Gergely Daróczi），介绍了信用违约模型，说明了如何使用copula对相关违约建模。

第8章“极值理论”（Zsolt Tulassay），给出了极值理论在保险和金融中的可能应用。你将学到如何对火灾损失分布的尾部拟合模型。然后用拟合模型计算在险价值（Value-at-Risk）和预期损失值。

第9章“金融网络”（Edina Berlinger，Gergely Daróczi），解释了金融网络在R中如何表示、模拟、可视化以及如何分析。我们将分析银行间借贷市场并学习如何系统化地检测重要的金融机构。

阅读本书之前的准备工作

本书提供的所有代码示例都应该在预装于计算机的R控制台上来运行。你可以免费下载软件并找到安装所有主要操作系统的指导。尽管本书没有包括高级主题，比如如何在整合发展环境中使用R，对于Emacs、Eclipse、vi或者Notepad++以及其他的编辑器有许多很好的插件，并且我们高度推荐你尝试RStudio，这是一款致力于R的免费开源的IDE。

除了R的安装版本，我们还会使用一些用户贡献的R包，这些包可以很容易地从CRAN（Comprehensive R Archive Network）进行安装。要安装一个R包，可以在R控制台使用install.packages命令，如下：

> install.packages('zoo')

安装之后，这个包需要在使用之前载入到当前的R会话中：

> library(zoo)

在R的主页，你能找到免费的入门文章和手册，但本书面向初学者，因此并不需要读者具备额外的R语言知识。

目标读者

本书为那些希望使用R来解决量化金融问题的读者而写。我们假定读者对金融有一定了解，但本书也会介绍金融理论。我们并不需要读者熟悉R，那些想开始学习R语言的读者会发现本书很有用，尽管我们没有给出完整的R语言概览，但说明了如何使用它的一部分来解决具体问题。即使你已经使用过R，也会惊讶于它所能应用的问题广度。

排版约定

在本书中，你会发现多种文本样式，用以区别不同种类的信息。这里举例说明其中一些类型，及其含义的解释。

文本中的代码字、数据库表格的名字、文件夹名称、文件名称、文件扩展名、路径、虚拟URL、用户输入以及推特的处理显示如下：“我们会运用一些forecast包中的方法。”

R代码块（通常是个函数体）安排如下。

logreturn <- function(x) {
    log(tail(x, -1) / head(x, -1))
}

当我们希望一个代码块的特定部分能吸引你的注意，相应的行或者项会设置为粗体。

logreturn tail(x, -1) / head(x, -1))
}

任何命令行的输入和输出的格式如下。

> pi
[1] 3.141593

其中R控制台中显示的“>”表示等待处理的命令。多行表达式会在第一行显示相同的符号，但其余诸行会在开头有一个“+”号，表示后面的R表达式尚待完成。

新术语和重要词汇以黑体表示。你在屏幕上看到的文字，例如，在菜单中或在对话框中，就像这样出现在文本中：“按下Next键，就翻到下一屏。”

警告或者重要的注解出现在这样的图标中。

提示或者技巧出现在这样的图标中。

读者反馈

我们始终欢迎读者的反馈。如果你对本书有任何想法，喜欢或者不喜欢什么，请让我们知道。读者的反馈对我们来说非常重要，这样我们才能出版读者最需要的图书。

一般性的反馈，请通过电子邮件发送到feedback@packtpub.com，请在邮件的主题中注明书名。

如果你精通某个领域并有兴趣写书或参与写书，请参考我们的作者指南www.packtpub.com/authors。

客户支持

现在，你已经是一位Packt图书的拥有者，我们会竭尽全力帮助你充分利用手中的书。

下载示例代码

可以使用你的账户从http://www.packtpub.com下载所有已购买的Packt图书的示例代码文件。如果你从其他地方购买了本书，可以访问http://www.packtpub.com/ support并注册，我们会通过电子邮件把文件发送给你。

勘误表

虽然我们已经竭力确保本书内容正确，但疏漏之处在所难免。如果你在我们的图书中发现错误，无论是文本还是代码，希望能通知我们，我们将不胜感激。这样做可以减少其他读者的困扰，帮助我们改进本书的后续版本。如果你发现任何错误，请访问http://www.packtpub.com/submit-errata，选择相应的图书，单击勘误表提交表单的链接，并输入详细的说明，然后提交。勘误一经核实，你的提交就被接受，此勘误将上传到本公司网站或者添加到现有勘误表。

要查看之前提交的勘误，登录https://www.packtpub.com/books/content/support并在搜索框中输入书名。请求的信息会在勘误部分出现。

侵权行为

互联网上的侵权材料是所有媒体都要面对的问题。在Packt，我们非常重视保护版权和许可证。如果你发现我们的作品在互联网上有任何形式的非法拷贝，请立即为我们提供网址或者网站名称，以便我们能够寻求补救。

请把可疑盗版材料发送到copyright@packtpub.com。

非常感谢你帮助我们保护作者以及宝贵的正版图书资源。

问题

如果你对本书内容有疑问，不管是哪个方面的，都可以通过questions@packtpub.com来联系我们，我们将尽最大努力解决。

第1章　时间序列分析

时间序列分析研究的是按时间顺序收集的数据。相邻的观测数据通常相互依赖。因此，时间序列分析的技术需要处理这种相依性。

本章的目标是通过一些特定应用来介绍一些常用建模技术。我们将看到如何使用R来解决现实中的这些问题。首先，我们考虑如何在R中存储和处理时间序列。接着，我们处理线性时间序列分析，并展现如何将它用于建模和预测房屋价格。其次，我们通过考虑长期趋势，使用协整的概念来改进基本的最小方差对冲比。最后，本章讲述如何将波动率模型运用于风险管理。

1.1　使用时间序列数据

用于存储时间序列数据的基本R类有vector、matrix、data.frame以及ts对象。但是，它们可以存储在这些对象中的数据类型相当有限。并且，这些表达方式提供的方法范围也很有限。不过幸运的是，同名的包中的特定对象，zoo、xts或timeSeries对象，对时间序列数据提供了更一般的表达形式。

对每个时间序列分析问题都创建时间序列对象是不必要的，但是复杂程度较高的分析则需要创建时间序列对象。你可以先将时间序列数据存储成向量形式，再计算数据的均值和方差，但如果你想用decompose对数据做季节分解，那就必须将数据存储在时间序列对象中。

下面的例子假定你使用了zoo对象，因为zoo对象是使用最广泛的包之一。在使用zoo对象之前，需要使用下面的命令安装并载入zoo包（如果你已经安装，那只需要载入它）。

>install.packages("zoo")
>library("zoo")

为了熟悉可用方法，我们使用存储在CSV文件aapl.csv中的苹果公司股票的日收盘价，创建了一个名为appl的zoo对象。表格的每一行包括一个日期和一个价格，两项通过逗号分隔。第一行包含了列名（Date和Close）。日期格式符合ISO8601推荐的基本标准符号（YYYY-MM-DD）。收盘价根据股票的拆分、股利以及相关改变进行调整。

小提示　下载示例代码
你对于在http://www.packtpub.com网站购买的所有Packt图书，都可以用自己在的账户从网站下载示例代码。如果你从其他途径购买了书籍，则可以访问http://www.packtpub.com/support并注册账号，示例代码会直接通过电子邮件发送给你。

使用下面的命令，可以从当前工作目录载入数据。

> aapl<-read.zoo("aapl.csv",
+　 sep=",", header = TRUE, format = "%Y-%m-%d")

为了初步了解数据，我们画出股票价格图形，并为整个图形设定一个标题（使用main参数）和对x轴和y轴标注了名称（分别使用xlab和ylab）。

> plot(aapl, main = "APPLE Closing Prices on NASDAQ",
+　 ylab = "Price (USD)", xlab = "Date")

使用下面的命令，我们可以提取时间序列开头部分或结尾部分。

> head(aapl)
2000-01-03 2000-01-04 2000-01-05 2000-01-06 2000-01-07 2000-01-10
　　 27.58　　　25.25　　　25.62　　　23.40　　　24.51　　　24.08
> tail(aapl)
2013-04-17 2013-04-18 2013-04-19 2013-04-22 2013-04-23 2013-04-24
　　 402.80　　 392.05　　 390.53　　 398.67　　 406.13　　 405.46

使用下面的命令，可以找出苹果股价在所有时间中的高点，和这个高点发生的日期。

> aapl[which.max(aapl)]
2012-09-19
　　694.86

当处理时间序列时，通常收益率更受关注，价格却不会。其原因是收益率通常平稳。因此我们会计算简单收益率或连续复合收益率（按百分比的形式）。

> ret_simple <- diff(aapl) / lag(aapl, k = -1) * 100
> ret_cont　 <- diff(log(aapl)) * 100

同时，我们也可以得到简单收益率的概括统计。在这里，我们使用coredata方法来表明我们仅仅关注股票价格，而非索引（日期）。

> summary(coredata(ret_simple))
　　 Min.　1st Qu.　Median　　Mean 3rd Qu.　　 Max.
-51.86000 -1.32500 0.07901 0.12530 1.55300 13.91000

可以看出，最大的单日损失是−51.86%。我们还可以使用下面的命令获得这个损失发生的日期。

> ret_simple[which.min(ret_simple)]
2000-09-29
 -51.85888

上网快速搜索可以发现，这个股价的剧烈变动缘于一个盈利预警的发布。我们可以画出直方图来加深理解日收益率的相关频率。对收益率数据进行分组时，我们可以使用break参数来指定分组数。

> hist(ret_simple, breaks=100, main = "Histogram of Simple Returns",
+　xlab="%")

我们也可以把分析限定于时间序列的一个子集（window）中。比如，苹果股价在2013年的最高点可以通过运行下面的命令的找到。

> aapl_2013 <- window(aapl, start = '2013-01-01', end = '2013-
12-31')
> aapl_2013[which.max(aapl_2013)]
2013-01-02
　　545.85

从风险管理的角度看，收益率分布的分位数很有意义。比如，我们使用简单的历史模拟法，可以很容易地确定置信水平为99%的一日在险价值（Value-at-Risk）。

> quantile(ret_simple, probs = 0.01)
　　　 1%
-7.042678

因此，在任意给定的一天中，收益率低于−7%的概率只有1%。但是如果这一天发生了这样的情形（每年大约会发生2.5次），7%将是最小的损失量。

线性时间序列的建模与预测

线性时间序列的一类重要模型是自回归单整移动平均（Autoregressive Integrated Moving Average，ARIMA）模型族，它由Box 和 Jenkins（1976年）提出。ARIMA模型假定了时间序列的当前值只依赖于自身的过去值和某些误差项的过去值。

根据Box和Jenkins的研究，建立ARIMA模型包含了以下3个阶段。

1．模型识别。

2．模型估计。

3．模型诊断检验。

模型识别的阶段包括了使用图方法或信息准则来确定试验模型的阶数（包含的过去值个数和过去误差项个数）。模型阶数确定之后需要估计模型参数，通常会使用最小二乘方法或者极大似然方法。最后，为了检查模型可能存在的缺陷，必须仔细检查拟合的模型。这个目的可以通过保证模型残差的行为符合白噪声的特点来实现，换句话说，残差不存在线性依赖。

1.2　对英国房屋价格建模并预测

除了zoo包，我们还会使用到forecast包的一些方法。如果你还没安装它，那就需要运行下面的命令来安装forecast包。

> install.packages("forecast")

接着我们运行下面的命令载入类。

> library("forecast")

首先，我们在时间序列对象zoo中存储月度房屋价格数据（来源：全英房屋抵押贷款协会）。

> hp <- read.zoo("UKHP.csv", sep = ",",
+　 header = TRUE,format = "%Y-%m", FUN = as.yearmon)

参数FUN对日期列调用给定的函数（as.yearmon，表示月度数据点）。为了确认指定as.yearmon真正存储了月度数据（每个周期12个子周期），我们可以查询数据序列的频率。

> frequency(hp)
[1] 12

结果表示，一个周期（称为年）中有12个子周期（称为月）。为了深入分析，我们再次计算数据的简单收益率。

> hp_ret <- diff(hp) / lag(hp, k = -1) * 100

1.2.1　模型识别和估计

我们使用forecast包提供的auto.arima函数，在一步中同时识别最优模型并估计参数。除了收益率序列（hp_ret），函数还使用了几个参数。通过指定stationary = TRUE，我们将搜索限定于非季节性模型。同样地，seasonal = FALSE将搜索限定于非平稳模型。此外，我们选择模型时，选择赤池信息量来度量模型的相对质量。

> mod <- auto.arima(hp_ret, stationary = TRUE, seasonal =FALSE,
+　 ic="aic")

为了确定拟合系数的值，我们可以查询模型的输出。

> mod
Series: hp_ret
ARIMA(2,0,0) with non-zero mean

Coefficients:
　　　　 ar1　　ar2 intercept
　　　0.2299 0.3491　　0.4345
s.e.　0.0573 0.0575　　0.1519

sigma^2 estimated as 1.105: log likelihood=-390.97
AIC=789.94　 AICc=790.1　 BIC=804.28

根据赤池信息量准则，看起来AR（2）模型拟合数据最好。为了视觉层次的确认，我们可以使用命令pacf画出偏自相关函数。图形显示了到滞后两阶的非零偏自相关，阶数为2的AR过程看来比较合适。AR（2）模型给出了两个AR系数、截距（如果模型包含AR项，截距实际上是均值）以及相应的标准误。在下面的例子中，因为水平为5%的置信区间没有包括0，所以这些统计量都在5%的水平上显著。

> confint(mod)
　　　　　　　2.5 %　　97.5 %
ar1　　　 0.1174881 0.3422486
ar2　　　 0.2364347 0.4617421
intercept 0.1368785 0.7321623

如果模型包含不显著的系数，我们可以使用带有fixed参数的arima函数重新估计模型，这相当于输入元素为0和NA的向量。NA表示相应的变量系数需要估计而0表示相应的变量系数需要设置为0。

1.2.2　模型诊断检查

一个快速验证模型的方法是运行下面的命令画出时间序列的诊断图。

> tsdiag(mod)

上述命令的输出在图1-1中显示。

图1-1　时间序列的诊断

标准化残差看来没有表现出波动率聚集，ACF图中的残差自相关并不显著，还有自相关的Ljung-Box检验的p值看起来很高，综上所以残差独立的原假设不能被拒绝，因此模型看来良好。

为了评估模型对样本数据的拟合良好程度，我们可以画出原始的月回报（细的黑色实线）与拟合值（粗的红色点线）的对比图形。

> plot(mod$x, lty = 1, main = "UK house prices: raw data vs.
fitted + values", ylab = "Return in percent", xlab = "Date")
> lines(fitted(mod), lty = 2,lwd = 2, col = "red")

输出显示在图1-2中。

图1-2　英国房屋数据

此外，我们可以计算精确性的常用度量。

> accuracy(mod)
ME　　　RMSE　　 MAE　　　MPE　　MAPE　　MASE
0.00120 1.0514　 0.8059　 -Inf　 Inf　　 0.792980241

这个命令返回平均误差、均方误差、平均绝对误差、平均百分比误差、平均绝对值百分比误差和平均绝对比例误差。

1.2.3　预测

为了预测接下来3个月的月收益率（2013年4～6月），我们使用下面的　　 命令。

> predict(mod, n.ahead=3)
$pred
　　　　　 Apr　　　 May　　　 Jun
2013 0.5490544 0.7367277 0.5439708
$se
　　　　　Apr　　　May　　　Jun
2013 1.051422 1.078842 1.158658

所以，我们预期在接下来的3个月中，平均房屋价格稍有增长，但标准误比较高，大约为1%。为了画出带有标准误的预测，我们可以使用下面的命令。

> plot(forecast(mod))

1.3　协整

协整的思想缘于Granger（1981年）提出的一个概念，后来由Engle和Granger（1987年）加以形式化。协整的思想是指寻找非平稳时间序列之间的一个线性组合，这个线性组合是一个平稳时间序列。因此，协整方法可以用于检测非平稳时间序列（比如价格）之间的稳定长期关系。

航空燃油的交叉对冲

航空公司很自然需要购买航空燃油。但由于航空燃油价格的波动很剧烈，大部分航空公司会将它们对航空燃油价格变化的风险敞口对冲掉一部分。如果市场中缺乏航空燃油OTC产品（译注：OTC产品指交易所场外柜台交易产品），航空公司会使用相关交易所交易的期货合约（比如，取暖油）来实现对冲。在下面的部分中，我们首先使用经典方法导出最优对冲比率，这种方法仅仅考虑两种价格之间短期波动。然后考察价格之间的长期稳定联系，进而改进最优对冲比。

首先，我们载入需要使用的包。urca包有一些有用的方法，可以用于单位根检验和估计协整关系。

> library("zoo")
> install.packages("urca")
> library("urca")

我们导入航空燃油和取暖油的月价格。

> prices <- read.zoo("JetFuelHedging.csv", sep = ",",
+　FUN =as.yearmon, format = "%Y-%m", header = TRUE)

现在我们仅仅考虑两种商品的短期行为（月价格改变）。通过拟合一个用取暖油价格变化来解释航空燃油价格的线性模型，可以推导出两种商品的最小方差对冲比率。线性模型中的beta系数就是最优对冲比。

> simple_mod <- lm(diff(prices$JetFuel) ~ diff(prices$HeatingOil)+0)

函数lm（用于线性模型）估计了航空燃油价格变动对取暖油价格变动的最佳拟合系数。+0项表示截距设置为0，意味着航空公司不持有现金。

> summary(simple_mod)
Call:
lm(formula = diff(prices$JetFuel) ~ diff(prices$HeatingOil) +
　　0)
Residuals:
　　 Min　　　 1Q　Median　　　3Q　　 Max
-0.52503 -0.02968 0.00131 0.03237 0.39602

Coefficients:
　　　　　　　　　　　 Estimate Std. Error t value Pr(>|t|)
diff(prices$HeatingOil) 0.89059　　0.03983　 22.36　 <2e-16***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.0846 on 189 degrees of freedom
Multiple R-squared: 0.7257, Adjusted R-squared: 0.7242
F-statistic: 499.9 on 1 and 189 DF, p-value: < 2.2e-16

结果得到对冲比为0.89059和残差标准误为0.0846。但是，交叉对冲并非完美无缺，最终得出的对冲组合结果仍然有风险。

现在，我们试图利用航空燃油和取暖油期货价格之间存在的长期关系来改进这一对冲比。我们使用下列命令画出两个价格序列（取暖油价格用红色），通过观察图形你可能已经猜出存在着这种长期关系。

> plot(prices$JetFuel, main = "Jet Fuel and Heating Oil Prices",
+　 xlab = "Date", ylab = "USD")
> lines(prices$HeatingOil, col = "red")

我们使用Engle和Granger的两步估计方法。首先，使用增强的Dickey-Fuller检验（augmented Dickey-Fuller test，ADF检验）对两个序列进行单位根检验（非平稳性）。

> jf_adf <- ur.df(prices$JetFuel, type = "drift")
> summary(jf_adf)
###############################################
# Augmented Dickey-Fuller Test Unit Root Test #
###############################################

Test regression drift
Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
　　 Min　　　 1Q　Median　　　3Q　　 Max
-1.06212 -0.05015 0.00566 0.07922 0.38086

Coefficients:
　　　　　 Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.03050　　0.02177　 1.401 0.16283
z.lag.1　　-0.01441　　0.01271　-1.134 0.25845
z.diff.lag　0.19471　　0.07250　 2.686 0.00789 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.159 on 186 degrees of freedom
Multiple R-squared: 0.04099, Adjusted R-squared: 0.03067
F-statistic: 3.975 on 2 and 186 DF, p-value: 0.0204

Value of test-statistic is: -1.1335 0.9865

Critical values for test statistics:
　　　1pct 5pct 10pct
tau2 -3.46 -2.88 -2.57
phi1　6.52　4.63　3.81

结果显示，因为检验统计量值−1.1335大于临界值−3.46，所以在1%的置信水平上不能拒绝非平稳（航空燃油时间序列包含一个单位根）的原假设。同样的结果对取暖油也成立（检验统计量是−1.041）。

> ho_adf <- ur.df(prices$HeatingOil, type = "drift")
> summary(ho_adf)

现在我们可以继续估计静态均衡模型，并使用ADF方法检验时间序列的残差是否平稳。请注意，目前的研究序列是上一步的估计结果，因此，我们现在必须使用不同的临界值［参见Engle和Yoo1987发表的论文］。

> mod_static <- summary(lm(prices$JetFuel ~ prices$HeatingOil))
> error <- residuals(mod_static)
> error_cadf <- ur.df(error, type = "none")
> summary(error_cadf)

得到的检验统计量是−8.912，而规模为200的样本在1%的置信水平上的临界值为−4.00。所以，我们拒绝非平稳的原假设。因此我们发现了两个协整变量并且可以进行第二步，两个协整变量指定一个误差修正模型（ECM）。ECM是一个动态模型，刻画了系统如何（以及多快）返回之前估计的静态均衡，这个静态均衡存储在变量mod_static中。

> djf <- diff(prices$JetFuel)
> dho <- diff(prices$HeatingOil)
> error_lag <- lag(error, k = -1)
> mod_ecm <- lm(djf ~ dho + error_lag)
> summary(mod_ecm)

Call:
lm(formula = djf ~ dho + error_lag + 0)

Residuals:
　　 Min　　　 1Q　Median　　　3Q　　 Max
-0.19158 -0.03246 0.00047 0.02288 0.45117

Coefficients:
　　　　　Estimate Std. Error t value Pr(>|t|)
Dho　 0.90020　　0.03238　27.798　 <2e-16 ***
error_lag -0.65540 0.06614　-9.909　 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.06875 on 188 degrees of freedom
Multiple R-squared: 0.8198, Adjusted R-squared: 0.8179
F-statistic: 427.6 on 2 and 188 DF, p-value: < 2.2e-16

通过考察航空燃油价格和取暖油价格之间存在的长期联系（协整），对冲比现在稍微提高了一点（0.90020），并且残差标准差显著地降低了（−0.06875）。误差项的系数为负（−0.65540）：两个价格之间原先较大的偏差有所修正，价格向它们的长期稳定关系移动。

1.4　波动率建模

正如我们之前所见，ARIMA模型常常用于过程的过去值已知时的条件期望建模。过去值已知的过程的条件方差是常数。真实世界的金融时间序列存在着波动性聚集和其他特点，换句话说，突发波动率打破了相对稳定的时期。

在这一节中，我们来考查GARCH时间序列模型。GARCH模型研究真实世界的（金融）时间序列的这个典型化事实——波动性聚集，并进一步运用这些模型预测在险价值（Value at Risk，VaR）。

1.4.1　风险管理的波动率预测

金融机构使用VaR来度量他们的活动风险，通常在10个工作日范围内计算置信水平为99%的临界值。这意味着在这10天内，预期实际损失只有1%的概率会超过VaR值。

我们载入zoo包并导入英特尔公司的月收益率数据，时间范围从1973年1月～2008年12月。

> library("zoo")
> intc <- read.zoo("intc.csv", header = TRUE,
+　 sep = ",", format = "%Y-%m", FUN = as.yearmon)

1.4.2　检验ARCH效应

收益率图形表明，在月收益率数据中可能存在ARCH效应。

> plot(intc, main = "Monthly returns of Intel Corporation",
+　 xlab = "Date", ylab = "Return in percent")

上面命令的输出在图1-3中显示。

图1-3　英特尔公司的月收益

我们可以使用统计假设检验来验证自己的想法。两种常用检验如下。

• 用于平方收益率（波动率的一种代理）自相关的Ljung-Box检验。
• Engle（1982年）提出的拉格朗日乘子（LM）检验。

首先，我们运行下面的命令，从而在平方收益率的前12阶滞后值上执行Ljung-Box检验。

> Box.test(coredata(intc^2), type = "Ljung-Box", lag = 12)

　 Box-Ljung test

data:　coredata(intc^2)
X-squared = 79.3451, df = 12, p-value = 5.502e-12

我们可以在1%的置信水平上拒绝原假设，原假设是平方收益率中不存在自相关。或者，我们可以使用FinTS包的LM检验，它输出相同的结果。

> install.packages("FinTS")
> library("FinTS")
> ArchTest(coredata(intc))

　 ARCH LM-test; Null hypothesis: no ARCH effects

data: coredata(intc)
Chi-squared = 59.3647, df = 12, p-value = 2.946e-08

两种检验都确定了英特尔的月收益率中存在ARCH效应。因此，收益率时间序列的建模应该使用ARCH或GARCH模型。

1.4.3　GARCH模型设定

GARCH（1,1）模型是GARCH模型中最常用的一种，也是一种最适于金融时间序列建模的模型。我们使用rugarch包提供的函数来设定模型、估计参数、回测以及预测。如果你还没有安装这个包，运行下面的命令。

> install.packages("rugarch")

然后，我们可以运行下面的命令来载入这个包。

> library("rugarch")

首先，我们需要使用函数ugarchspec设定模型。对于一个GARCH（1,1）模型，我们需要设置garchOrder为c(1,1)。而且均值模型（mean.model）是一个白噪声过程，因此等同于armaOrder = c(0,0)。

> intc_garch11_spec <- ugarchspec(variance.model = list(
+　 garchOrder = c(1, 1)),
+　mean.model = list(armaOrder = c(0, 0)))

1.4.4　GARCH模型估计

利用极大似然方法的系数的精确拟合，由使用了模型指定和收益率数据作为输入的ugarchfit函数来完成。

> intc_garch11_fit <- ugarchfit(spec = intc_garch11_spec,
+　data = intc)

其他参数的使用请参见ugarchfit的帮助文档。拟合模型的输出（使用命令intc_garch11_fit）展示了有用信息，比如最优参数值、对数似然函数值，以及信息准则。

1.4.5　回测风险模型

检验模型表现的有效方法是历史回测。在回测风险模型时，我们对比整个时期的真实收益率和估计的VaR。如果收益率比VaR损失更大，我们得到一次VaR突破（VaR exceedance）。在我们的例子中，一次VaR突破应该仅仅发生在1%的情形中（因为我们设定了置信水平为99%）。

函数ugarchroll对一个设定的GARCH模型（在这里，这个模型是intc_garch11_spec）执行历史回测。我们指定回测如下。

• 使用的收益率数据存储在zoo对象intc中。
• 回测的起始点（n.start）应该是序列开始后的第120个月（即1983年1月）。
• 每月都需要重新估计模型（refit.every = 1）。
• 我们使用移动（moving）窗口来估计。
• 我们使用一个混合（hybrid）的解决方法。
• 我们希望计算VAR在99%的尾部置信水平上（VaR.alpha = 0.01）的临界值（calculate.VaR = TRUE）。
• 我们希望保留估计的系数（keep.coef = TRUE）。

下面的命令显示了满足上述所有要求的回测。

> intc_garch11_roll <- ugarchroll(intc_garch11_spec, intc,
+ n.start = 120, refit.every = 1, refit.window = "moving",
+ solver = "hybrid", calculate.VaR = TRUE, VaR.alpha = 0.01,
+ keep.coef = TRUE)

我们可以使用report函数检查回测报告。通过把这个参数的type参数设定为VaR，这个函数对突破值执行无条件和有条件覆盖检验。VaR.alpha是尾部概率，conf.level是置信水平，条件覆盖的假设检验基于此建立。

> report(intc_garch11_roll, type = "VaR", VaR.alpha = 0.01,
+conf.level = 0.99)
VaR Backtest Report
===========================================
Model:　　　　　　sGARCH-norm
Backtest Length:　 312
Data:

==========================================
alpha:　　　　　　　1%
Expected Exceed:　　 3.1
Actual VaR Exceed:　　5
Actual %:　　　　　 1.6%

Unconditional Coverage (Kupiec)
Null-Hypothesis:　 Correct Exceedances
LR.uc Statistic:　 0.968
LR.uc Critical:　　　6.635
LR.uc p-value:　　　0.325
Reject Null:　　　NO

Conditional Coverage (Christoffersen)
Null-Hypothesis:　 Correct Exceedances and
　　　　　　　 Independence of Failures
LR.cc Statistic:　 1.131
LR.cc Critical:　　　9.21
LR.cc p-value:　　　0.568
Reject Null:　　　 O

Kupiec的无条件覆盖方法比较了给定VaR尾部概率时，预期突破值数目和实际突破值数目，而Christoffersen检验方法则是一种对无条件覆盖和突破值的独立性的联合检验。在我们的例子中，尽管预期突破有3次但实际发生了5次，我们不能拒绝突破是正确并且独立的原假设。

回测表现的图形也很容易生成。首先，使用ugarchroll对象提取的预测VaR创建一个zoo对象。

> intc_VaR <- zoo(intc_garch11_roll@forecast$VaR[, 1])

我们仍然使用这个“zoo”对象，通过rownames（年和月）重写这个对象的index属性。

> index(intc_VaR) <- as.yearmon(rownames(intc_garch11_roll@forecast$VaR))

对同时存储在ugarchroll对象中的真实收益率，我们加以同样的处理。

> intc_actual <- zoo(intc_garch11_roll@forecast$VaR[, 2])
> index(intc_actual) <-
as.yearmon(rownames(intc_garch11_roll@forecast$VaR))

现在，我们可以使用下面的命令，画出VaR对比英特尔真实收益率的图形。

> plot(intc_actual, type = "b", main = "99% 1 Month VaR Backtesting",
+　 xlab = "Date", ylab = "Return/VaR in percent")
> lines(intc_VaR, col = "red")
> legend("topright", inset=.05, c("Intel return","VaR"), col =
c("black","red"), lty = c(1,1))

图1-4中显示了上述命令行的输出。

图1-4　99% 1月VaR回测

1.4.6　预测

我们现在有理由相信风险模型运行正常，我们也可以生成VaR预测。函数ugarchforecast选取以下两个参数，一个是拟合的GARCH模型（intc_garch11_fit），另一个是应该产生预测的周期数（n.ahead = 12，即12个月）。

> intc_garch11_fcst <- ugarchforecast(intc_garch11_fit, n.ahead = 12)

我们可以通过查询以下命令行显示的预测对象，来预期未来结果。

> intc_garch11_fcst
*------------------------------------*
*　　　　　GARCH Model Forecast　　　　*
*------------------------------------*
Model: sGARCH
Horizon: 12
Roll Steps: 0
Out of Sample: 0

0-roll forecast [T0=Dec 2008]:
　　　Series Sigma
T+1　0.01911 0.1168
T+2　0.01911 0.1172
T+3　0.01911 0.1177
T+4　0.01911 0.1181
T+5　0.01911 0.1184
T+6　0.01911 0.1188
T+7　0.01911 0.1191
T+8　0.01911 0.1194
T+9　0.01911 0.1197
T+10 0.01911 0.1200
T+11 0.01911 0.1202
T+12 0.01911 0.1204

波动率（sigma）的一步预期是0.1168。因为我们假定了正态分布，置信水平为99%的VaR可以使用标准正态分布的99%分位数（输入qnorm(0.99)）来计算。因此对下一个周期，一个月的99%VaR就是qnorm(0.99)*0.1168 =0.2717。结果，月收益率高于−27%的概率是99%。

1.5　小结

在本章中，我们将R运用于时间序列分析中一些精选出来的问题。本章涵盖了时间序列数据的不同表达方法，使用ARMA模型预测房屋价格，使用协整关系改进基本的最小方差对冲比，以及使用GARCH模型进行风险管理。在下一章，你将学习如何使用R来构建最优投资组合。

第2章　投资组合优化

到现在为止，我们已经熟悉了R语言的基础。我们知道如何去分析数据、调用它的内置函数并把它们运用到我们选择的时间序列分析问题上。在本章中，我们既运用并扩展这些知识来讨论一种重要的实践应用，即投资组合分析，换句话说也叫证券选择。这一节涵盖了投资组合优化背后的思想：数学模型和理论求解。为了提高编程技巧，我们使用真实数据解决一个现实中的问题，并逐行实施算法。同时，我们也在相同的数据集上使用预先写好的R包。

想象我们生活在一个热带岛屿，只有100美元可以投资。这个岛上的投资机会相当有限。我们可以把全部资金投资到冰淇淋上或雨伞上。收益取决于天气，如表2-1所示。

表2-1　　天气与收益数据

假定天气是晴天或是雨天的概率相同。如果我们不能预知天气或者改变天气，这两种选择的概率显然相等，我们投资于其中任何一种，都会得到5%的预期收益率［(0.5×120+0.5×90)/100−1=0.05］。

如果我们可以把资金分配在冰淇淋和雨伞之间，那该如何划分资金？我们应该在两种备选中各投资50美元。无论会发生什么，我们会在一种资产上赚到45美元而在另一种资产上赚到60美元，因此，这个组合没有风险。预期收益率仍然是5%，但因为(45+60)/100−1=0.05，现在的收益获得了保障。

这个例子抓住了投资组合优化的主要概念（正是因为这个理论，哈里•马科维茨在1990年获得诺贝尔奖）。通过考虑投资产品之间的相关性，我们可以在保持想要的预期收益率不变的前提下减少投资组合的风险（在这个例子里由方差来表示）。

为了得到精确的数学表达式，令X和Y是随机变量，各自方差有限并分别为和。它们的凸组合或仿射组合的方差显示在下面的二次函数中。

对于它们相关系数的不同值，这个二次函数如图2-1所示。

图2-1　投资组合的方差

当且仅当X和Y的相关系数为−1，并且X和Y的方差相同时，方差（作为风险的一种度量）可以完全消除掉（译者注，原文此处有误）。否则，我们会在后面的“定理（拉格朗日）”这一节中看到，权重最优的投资组合的方差取决于（以一种完全不平凡的方式）所有3个参数，分别是、以及。

2.1　均方差模型

马科维茨[Markowitz,H.M.（March 1952）]提出的均方差模型，其实是冰　 淇淋/雨伞交易模型在高维上的体现。为了得到数学公式，我们需要一些相关的定义。

定义的解释如下。

• 资产X_i意味着一个方差有限的随机变量。
• 投资组合意味着资产的组合：，其中并且 。这个组合可以是仿射的或者凸的。在凸组合的情况下，所有的资产权重都非负。
• 最优化意味着选择最佳的w_i系数（权重）的过程，以便使投资组合满足我们的需要（也就是说，在给定的预期收益率水平上选择最小的风险，或者在给定的风险水平上选择最高的预期收益率，等等）。

令X₁, X₂, …，X_n是方差有限的随机收益率变量，是它们的协方差矩阵，并且。

我们会关注下列的最优化问题。

（1）

（2）

（3）

（4）

（5）

很明显，w^TQw是投资组合的方差，wr是预期收益率。对于权重和，我们有，这意味着我们愿意投资一个单位的资金。（我们也可以考虑增加的条件，这意味着不允许做空。）下面各点详细地解释了问题。

• 第一个问题是寻找风险最小的投资组合。如果考虑无风险资产，将资产全部配置在无风险资产上，最小风险就是零了，非平凡就是为了避免这样的情况。
• 第二个问题是在给定的方差水平上选择最大的预期收益率。
• 一个稍微不同的方法是在一个想要的预期收益率水平上寻找方差最小的投资组合。
• 第四个问题是最大化一个简单的效用函数，return−λ×risk，其中λ是风险容忍的系数，它以一个任意数字表示我们对风险的态度。这个问题实际上等同于第一个问题。
• 在第五个问题中，Y是第n+1个资产（比如说，一个指数），是一种我们不能买或者不想买，但需要复制的资产。其他类似的问题可以按同样的方式公式化。

很明显，第二个问题是一个带有二次约束的线性优化问题，而其他所有问题是带有线性约束的二次函数最优化问题。随后我们可以看到，这是一个相当重要的区别，因为线性约束容易处理而二次约束很难处理。在接下来的两节中，我们会专注于这些问题的复杂性和可能的解。

2.2　解的概念

在过去的50年中，人们开发了许多表现良好的伟大算法用于数值优化，尤其用于二次函数的情形。正如我们在上一节中看到的，我们只有二次函数和约束。所以这些方法（它们也在R中实现）可以用于最糟的情形（如果没有任何更好的情形）。

然而，数值优化的细节讨论超出了本书范围。幸运的是，在线性约束和二次函数约束的特殊情况下，这些方法并不必要。我们可以使用18世纪提出的拉格朗日定理。

（拉格朗日）定理

如果f:Rⁿ→R和g=(g₁, g₂, …，g_m):Rⁿ→R^m（其中m<n）有连续的偏导数，并且是使得g(x)=0的f(x)的相对极值点，其中限定rank(g'(a))=m。

那么，这里存在系数λ₁, λ₂, …, λ_m使得。换句话说，这个函数所有的偏导数都是0（Bertsekas Dimitri P. (1999))。

在我们的例子里，这个条件也是充分的。二次函数的偏导数是线性的，所以最优化产生了求解一个线性方程组的问题，（不像数值优化）这是一个高中水平的任务。

我们来看一下，如何运用它来解决第三个问题。

可以表明，这个问题等价于下面的线性方程系统。

（两行和两列加入到协方差矩阵中，所以我们有条件同时决定两个拉格朗日乘子。）我们可以期望这个系统有唯一解。

需要强调的是，我们利用拉格朗日定理得到的结果已经不再是一个最优化问题。就像是在一维中，最小化二次函数需要求导和解线性方程组，但从复杂性的角度看，求解并不复杂。现在，我们来看看如何处理收益率最大化的问题。

容易看出，拉格朗日函数对于λ的导数是约束本身。

为了看出这一点，我们取L的导数。

由此产生了一个非线性方程，它的求解与其说是科学不如说是艺术。

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是一种公共服务，提供了轻松获取大批定量数据的途径。这个超链接［（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

接着，我们来建立拉格朗日定理指定的线性等式系统的左侧：，其中我们按行（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（在最小方差函数的参数中指定）。这是线性系统右侧向量的最后一个值。

> 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中。

图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所示的一样。

图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中。

图2-4　有效前沿

2.4　切线组合和资本市场线

当组合中加入一个无风险资产R，会发生什么？如果σR=0并且X是任意的一个风险投资组合，那么Var(α_R+(1−α)X)=(1−α)²Var(X)，并且显然也有E(αR+(1−α)X)= αE(R)+(1−α)E(X)。这意味着这些组合在均值-标准差平面上形成了一条直线。位于这条直线上的任何投资组合都可以通过投资于R和X来得到。很明显，X的最佳选择位于这条直线与有效边界的切点。这个切点叫作市场组合或者切点组合，而风险资产的有效前沿在这个点的切线叫作资本市场线（Capital Market Line，CML），它包含了该情形下所有资产的有效投资组合。我们解决的最后一个关于均值－方差模型的问题是，如何决定市场投资组合（或等价的，CML）。

通过修正方差最小化的代码，我们可以轻松地解决这个问题。首先，如果我们想加入无风险资产，就在协方差矩阵中加入了全部为0的一行和一列（其中n是总资产个数，包括无风险资产在内）。

> n <- 6; mu <- 0.005
> Q <- cbind(cov(return), rep(0, n - 1))
> Q <- rbind(Q, rep(0, n))

并在收益率向量中加入无风险资产。

> r <- c(colMeans(return), rf)

然后，我们可以使用新的协方差矩阵和新的收益率向量来决定最优投资组合的权重，接着基于“使用真实数据工作”这一节描述的minvariance代码来移除第n个资产。

> Q <- rbind(Q, rep(1, n), r)
> Q <- cbind(Q, rbind(t(tail(Q, 2)), matrix(0, 2, 2)))
> b <- c(rep(0, n), 1, mu)

得到下面的中间结果。

> round(Q, 6)
　　　　 AAPL　　GOOG　　　MSFT　　 IBM　　　　　　　T　　　　　　　　r
AAPL　0.000630 0.000338　0.000249 0.000233　0.000218 0e+00 1　0.000748
GOOG　0.000338 0.000462　0.000226 0.000186　0.000182 0e+00 1　0.000008
MSFT　0.000249 0.000226　0.000341 0.000178　0.000177 0e+00 1 -0.000236
IBM　 0.000233 0.000186　0.000178 0.000240　0.000157 0e+00 1　0.000439
T　　 0.000218 0.000182　0.000177 0.000157　0.000283 0e+00 1 -0.000179
　　　0.000000 0.000000　0.000000 0.000000　0.000000 0e+00 1　0.000100
　　　1.000000 1.000000　1.000000 1.000000　1.000000 1e+00 0　0.000000
r　　 0.000748 0.000008 -0.000236 0.000439 -0.000179 1e-04 0　0.000000
> b
[1] 0.000 0.000 0.000 0.000 0.000 0.000 1.000 0.005

解方程后，市场组合的结果如下。

> w <- solve(Q, b)
> w <- head(w, -3)
> w / sum(w)
　　　AAPL　　 GOOG　　　MSFT　　　　IBM　　　　 T
-10.154891 4.990912 12.347784 -18.010579 11.826774

2.5　协方差矩阵中的噪声

当我们优化投资组合时，其实我们并没有真实的协方差矩阵和预期收益率向量（它们是均值-方差模型的输入量）。我们使用观测来估计它们，因此Q、r以及模型的输出都是随机变量。

如果不深入细节，我们可以说模型中会产生惊人的巨大不确定性。尽管有强大数定律的保证，最优的投资组合权重会不时地在±200%之间变动。幸运的是，如果我们掌握有几年的数据（日收益率），测量风险的相对误差就仅为20%～25%。

2.6　如果方差不够用

方差可以很容易地度量风险，但也存在一些缺陷。例如，在使用方差时，收益率中的正向变化也会被视为风险的增加。因此，人们开发了一些更复杂的风险度量方法。

比如下面的这个简例，主要关于多种方法运用于之前描述过的IT_return资产，它的目的是快速回顾fPortfolio包提供的选择。

> Spec <- portfolioSpec()
> setSolver(Spec) <- "solveRshortExact"
> setTargetReturn(Spec) <- mean(colMeans(IT_return))
> efficientPortfolio(IT_return, Spec, 'Short')
> minvariancePortfolio(IT_return, Spec, 'Short')
> minriskPortfolio(IT_return, Spec)
> maxreturnPortfolio(IT_return, Spec)

这些R表达式返回了不同的组合权重，它们的计算方法很多，但没有在这个介绍性的章节中讨论。细节请参考这个包的捆绑文档，比如?fportfolio，以及列在参考文献中的相关文章和图书的章节。

2.7　小结

本章的主题是投资组合优化。我们先介绍了主要的思想，接着介绍了马可维茨模型和它的数学公式。我们讨论了可能的求解方法，并运用一个简单的算法来展示这些方法如何运用于真实数据。我们还使用预先写好的R包来解决相同的问题。我们广泛讨论了其他重要的主题，比如市场组合、协方差矩阵在估计中的不确定性，以及超越方差的风险度量。我们希望它们是这些主题的有益开始，能够鼓励你进行深入研究或者接着阅读下一章，它讲述一个相关的主题——资产定价模型。