详情
这是一本介绍量化金融的书，书中运用R语言这一工具教会读者实现定量分析，在实际的交易环境中实现更好的操作。书中不仅涉及很多与量化金融相关的金融学专题，同时也详细介绍了R语言相关的编程技巧。本书适合一切有志于通过计算机编程更好地优化交易成果的人士。
书名：量化金融R语言高级教程
ISBN：978-7-115-44982-5
本书由人民邮电出版社发行数字版。版权所有，侵权必究。
您购买的人民邮电出版社电子书仅供您个人使用，未经授权，不得以任何方式复制和传播本书内容。
我们愿意相信读者具有这样的良知和觉悟，与我们共同保护知识产权。
如果购买者有侵权行为，我们可能对该用户实施包括但不限于关闭该帐号等维权措施，并可能追究法律责任。
• 著 [匈牙利] Edina Berlinger 等
译 高 蓉
责任编辑 胡俊英
• 人民邮电出版社出版发行 北京市丰台区成寿寺路11号
邮编 100164 电子邮件 315@ptpress.com.cn
• 读者服务热线：(010)81055410
反盗版热线：(010)81055315
Copyright ©2015 Packt Publishing. First published in the English language under the title Mastering R for Quantitative Finance.
All rights reserved.
本书由英国Packt Publishing公司授权人民邮电出版社出版。未经出版者书面许可，对本书的任何部分不得以任何方式或任何手段复制和传播。
版权所有，侵权必究。
R语言是用于统计分析、绘图的语言和操作环境，是属于GNU系统的一个自由、免费、源代码开放的软件。它是一个用于统计计算和统计制图的优秀工具。
本书通过13章的内容向读者详细介绍了使用R语言实现量化金融的方方面面。本书包括实证金融（第1～4章）、金融工程（第5～7章）、交易策略优化（第8～10章）和银行管理（第10～13章）等主题。
本书的目标读者是那些既熟悉基本金融概念又具有一定编程能力的人。通过阅读本书，读者可以了解R语言与量化金融相关的各类知识和编程技巧。
自20世纪50年代马可维茨开创投资组合理论以来，金融学开始与数学模型和数据分析紧密融合，逐渐形成了高度数理化的现代金融学。时至今日，无论是对传统的金融学理论，还是对如日中天的量化投资，又或是对方兴未艾的金融科技（Fintech），数学模型与数据科学源源不断地为这些金融理论和实践提供着发展工具，一起构建了复杂深奥、丰富多彩的量化金融世界。由于数学模型与数据分析本身的复杂性，量化金融技术一直被公众视为数学家或华尔街的“火箭科学家”的秘密武器。
近年来，开源软件得到了迅速推广，促进了各行各业对数据科学的学习与应用。R就是这样一种开源数据分析工具。它灵活易用，数据分析能力强大，从技术上降低了量化金融技术的学习和应用难度，目前在世界范围内被广泛地采用，在中国也拥有大量用户。但是，可以结合量化金融与数据科学的图书，在市场上依然很少，而本套书正是基于这种需求而作的关于量化金融的学习教程。本套书选择R语言作为工具，通过R语言学习量化金融，帮助读者快速掌握知识并加以实践应用，为各种金融问题提供实践解决方案，同时还可使用第三方贡献的免费R包。
本套书内容广泛、取舍精当，涵盖了实证金融、金融工程、交易策略与银行管理等内容，可以帮助读者全面学习金融理论与技术。本套书既包含金融时间序列分析、资产定价和期权定价等理论知识，也包括投资组合管理、信用风险管理等实践知识，还涉及金融学的前沿领域以及2008年金融危机后发展起来的金融网络分析理论。本套书包括两本，分别是《量化金融R语言初级教程》和《量化金融R语言高级教程》。学习《量化金融R语言初级教程》不需要读者精通R和金融理论，但需要对金融领域具有一定了解。学习《量化金融R语言高级教程》要求读者具备一定的R编程能力和金融学基本概念，建议读者在学习“初级教程”的基础上再学习“高级教程”。
作为一名金融学教师，我认为本套书是通向量化金融世界的一把“金钥匙”，可以有效地促进金融知识和计算技术的传播。我很荣幸能够承担本书的翻译工作。在此特别感谢胡俊英编辑，胡编辑认真专业的审核工作，有力保证了本书翻译的如期完成，我亦受益匪浅。同时特别感谢天津理工大学的李茂老师，书中多处的专业词汇敲定，都来自与李茂老师的讨论结果。当然，文责自负。同时，感谢杭州电子科技大学2016年高等教育研究资助项目YB201631“投资学教学与R软件应用”的支持。
——高蓉，2017年3月
高蓉，任教于杭州电子科技大学经济学院。博士毕业于南开大学经济学院，本科毕业于南开大学数学科学学院。研究领域包括资产定价、实证金融、数据科学应用。已出版教材《实验投资学》，译著《数据科学入门》，发表学术论文数篇。
本书是我们的前一本书《量化金融R语言初级教程》（Introduction to R for Quantitative Finance）的续作。本书是为那些希望学习R语言来建立更高级量化金融模型的读者而写的。本书包括实证金融（第1～4章）、金融工程（第5～7章）、交易策略优化（第8～10章）和银行管理（第10～13章）等主题。
第1章，“时间序列分析”讨论了一些重要的概念，比如，协整（结构化）、向量自回归模型、脉冲-响应函数、基于非对称GARCH模型的波动率建模以及信息影响曲线。
第2章，“因素模型”讲述了如何建立多因子模型并实施。借助主成分分析方法，可以识别出5个独立因子解释资产回报率。为了给出例证，我们基于真实市场数据重新建立了Fama-French模型。
第3章，“成交量预测”讲述了日内成交量预测模型，及其基于DJIA指数数据的R语言实现。模型没有使用成交率，而是使用了换手率，并从动态成分中分解出季节成分，进而分别预测了这两个成分。
第4章，“大数据—高级分析”使用R访问开源数据，并在大数据上执行多种分析。为了给出例证，我们将K-均值聚类和线性回归运用到了大数据上。
第5章，“FX衍生品”推广了用于衍生品定价的Black-Scholes模型。Margrabe公式是对Black-Scholes 模型的一种扩展，可以用来编程实现对股票期权、货币期权、外汇期权和quanto期权的定价。
第6章，“利率衍生品和模型”给出了利率模型和利率衍生品的概览。使用Black模型对上限和上限单元定价，并给出了利率模型，如Vasicek模型和CIR模型。
第7章，“奇异期权”介绍了奇异期权，解释了它们与普通香草期权的关系，并对任意衍生品定价函数估计了它们对应的希腊字母表示，并且更详细地考察了一种特殊的奇异期权——无触发二元期权。
第8章，“最优对冲”分析了衍生品对冲过程中的一些应用问题，这些问题由组合在离散时间上重新安排和交易成本引起。为了找到最优对冲策略，使用了不同的数值算法。
第9章，“基本面分析”研究了如何基于基本面建立交易策略。为了选出收益最佳的股票，一方面，根据公司过去的表现创建它们的聚类；另一方面，借助决策树区分出表现卓越的股票。基于这些结果，定义股票选择的规则并实施了回测。
第10章，“技术分析、神经网络和对数优化组合”给出了技术分析和一些相应策略的概述，如神经网络和对数优化组合；并以动态的设定研究了预测单资产（比特币）价格、优化交易择时以及（NYSE股票）组合配置的问题。
第11章，“资产和负债管理”讲述了如何使用R语言支持银行的资产和负债管理。关注点在于数据生成、利率风险的度量和报告、流动性的风险管理以及无到期存款行为的建模。
第12章，“资本充足率”总结了巴塞尔协议的原则，而且借助历史的、delta-正态的、蒙特卡洛模拟的方法计算了在险价值（VaR），用以确定银行的资本充足性。同时也对信用和操作风险的特定问题有所涉猎。
第13章，“系统风险”讲述了两种方法，这些方法基于网络理论，并有助于识别系统性重要的金融机构，分别是核心-边缘模型和传染模型。
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('Quantmod')这个包安装完毕之后，需要先装载在当前的R会话中，再进行使用：
> library (Quantmod)在R的主页，你可以找到免费的入门文章和手册。
本书的目标读者是那些既熟悉基本金融概念又具有一定编程能力的人。但是，即使你具备量化金融的基础，或者你已经拥有一些R的编程经验，本书同样可以给你新的启示。即使你已经是相关领域的专家，本书也能帮助你迅速切入其他领域。不过，如果你希望完美地掌握各章的内容，需要具备中级水平的量化金融知识，并且还需要掌握关于R的相关知识。以上这些技能可以从本书的前作——《量化金融R语言初级教程》中获得。
在本书中，你会发现一些不同的文本样式，用以区别不同种类的信息。这里举例说明其中一些样式，以及它们的含义。
命令行的输入和输出的格式如下：
#generate the two time series of length 1000
set.seed(20140623) #fix the random seed
N <- 1000 #define length of simulation
x <- cumsum(rnorm(N)) #simulate a normal random walk
gamma <- 0.7 #set an initial parameter value
y <- gamma * x + rnorm(N) #simulate the cointegrating series
plot(x, type='l') #plot the two series
lines(y,col="red")新术语和重要的词以黑体表示。你在屏幕上看到的文字，例如，在菜单或对话框中，就像这样出现在文本中：“另一种有用的可视化练习是看对数尺度上的密度。”
警告或者重要的注解出现在这样的图标中。
提示或者技巧像这样出现。
我们始终欢迎读者的反馈。如果你对本书有任何想法，喜欢或者不喜欢什么，请让我们知道。读者的反馈对我们来说非常重要，这样我们才能出版读者最需要的图书。
一般性的反馈，请通过电子邮件发送到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来联系我们，我们将尽最大努力来解决问题。
Edina Berlinger拥有布达佩斯考文纽斯大学的博士学位。她是一位助理教授，讲授公司财务、金融学和金融风险管理。她还是大学金融系的领导，也是匈牙利科学院金融分会的主席。她的专业领域涉及信贷系统、风险管理以及最近的网络分析。她已经领导过几个研究项目：学生贷款设计、流动性管理、异质代理模型和系统风险。
“本工作由匈牙利科学院的动量项目（LP-004/2010）支持。”
Ferenc Illés拥有罗兰大学的数学硕士学位。毕业之后的一些年中，他开始研究精算和金融数学，而且他即将开始在布达佩斯考文纽斯大学的博士学习。最近几年，他在银行业工作。目前他正在开发使用R的统计模型。他的兴趣与大型网络以及计算复杂性有关。
Milán Badics拥有布达佩斯考文纽斯大学的硕士学位。现在，他是一名博士生，并且是PADS博士奖学金项目的成员。他讲授金融计量学，而且他的研究主题是使用数据挖掘方法的时间序列预测、金融信号处理以及利率模型的数值敏感分析。2014年5月，他在由匈牙利证券交易所组织的X. Kochmeister奖项的竞赛中获胜。
Ádám Banai从布达佩斯考文纽斯大学得到投资分析和风险管理的硕士学位。他加入了匈牙利国家银行（Magyar Nemzeti Bank，MNB，匈牙利的中央银行）的金融稳定性部门。从2013年起，他成为金融系统分析理事会（MNB）应用研究和压力测试部门的领导。自2011年起，他也是布达佩斯考文纽斯大学的博士生。他的主要研究领域是偿付能力压力测试、资金流动性风险和系统风险。
Gergely Daróczi是一位狂热的R包开发者，并且是一家位于Rapporter的R网络应用公司的创始人和CTO。他同时也在攻读社会学博士学位，并且目前作为R开发者领导在洛杉矶的CARD.com工作。如果算上讲授统计学和从事数据分析项目的几年时间，他大约已经有10年的R编程环境的工作经验。Gergely是《量化金融R语言初级教程》（Introduction to R for Quantitative Finance）的合著者，目前除了一些关于社会科学的杂志文章和报告，他同时还忙于另一本Packt出版社的图书《精通R语言数据分析》（Mastering Data Analysis with R）。他对那本书的贡献是审阅并负责R源代码的格式。
Barbara Dömötör是布达佩斯考文纽斯大学金融系的一名助理教授。在2008年开始博士学习之前，她曾为多家跨国银行工作。她的博士论文与公司的套期保值有关。她撰写了关于公司财务、金融风险管理和投资分析的讲义。她的主要研究领域是公司财务、金融风险管理和公司的套期保值。
Gergely Gabler自2014年起是匈牙利国家银行（MNB）金融监管单位的商业模型分析部门领导。在这之前，自2008年起，他曾经是匈牙利Erste银行宏观经济研究部门的领导人。他在2009年毕业于布达佩斯考文纽斯大学，并获得金融数学的硕士学位。自2010年起，他在布达佩斯考文纽斯大学任客座讲师，同时也在MCC学院做高等研究的讲座。他预计会在2015年结束CFA考试，并成为一名持证人。
Dániel Havran是一名匈牙利科学院经济和区域研究中心经济研究所的博士后研究人员。他同时作为布达佩斯考文纽斯大学兼职助理教授，在那里他讲授公司财务（本科、博士）以及信用风险管理（硕士）。在2011年，他获得了布达佩斯考文纽斯大学的经济学博士学位。
“我非常感谢匈牙利科学院博士后奖学金计划的支持。”
Péter Juhász拥有布达佩斯考文纽斯大学的工商管理博士学位，同时也持有CFA证书。作为一名助理教授，他讲授公司财务、商业估值、Excel的VBA编程以及沟通技巧。他的研究领域涉及无形资产的估值、商业表现分析和建模以及政府采购和体育管理。他曾写过一些文章和书的某些章节，主要关于匈牙利公司的财务表现。同时，他也定期为中小企业服务，而且在安永商业学院的EMEA（欧洲、中东和非洲）区域任高级培训师。
István Margitai是CEE（中东欧）区域一家主要银行集团的资产负债管理团队的分析师。他主要处理方法论问题、产品建模以及内部转移定价等主题。在2009年，他开启了在匈牙利的资产负债管理的职业生涯，并收获了战略流动性管理和流动性计划的经验。他在布达佩斯考文纽斯大学主修投资和风险管理。他的研究兴趣是银行业的微观经济学、市场微观结构以及订单驱动市场的流动性。
Balázs Márkus从事金融衍生品工作已经超过10年。他曾经交易过多种类型的衍生品，从碳互换到国债期货的期权。他是布达佩斯Raiffesien银行外汇衍生品部门的领导。他是Pallas Athéné Domus环境科学基金会顾问委员会的一员、匈牙利国家银行的兼职分析师以及一家小型的证券自营和顾问公司Nitokris有限公司的常务董事。目前，他正在布达佩斯考尔纽斯大学攻读动态对冲作用的博士学位，同时他还是那里的一名教学助理。
Péter Medvegyev拥有布达佩斯Marx Károly大学的经济学硕士学位。在1977年毕业之后，他开始了匈牙利管理发展中心的顾问工作。他在1985年获得了经济学博士学位。自1993年开始，他为布达佩斯考文纽斯大学数学系工作。他在考文纽斯大学的教学经历涵盖随机过程、数理金融以及其他多门数学专业课。
Julia Molnár是布达佩斯考文纽斯大学的一名博士学位候选人。她的主要研究兴趣包括金融网络、系统风险以及零售银行业的金融技术创新。自2011年起，她为McKinsey & Company工作，在那里她参与了银行业领域的多项数字和创新研究。
Balázs Árpád Szűcs是布达佩斯考文纽斯大学的金融学博士生，并同时在该大学的金融系任研究助理。他拥有投资分析和风险管理的硕士学位。他的研究兴趣包括最优执行、市场微观结构和日内交易量预测。
Ágnes Tuza拥有布达佩斯考文纽斯大学的应用经济学学位，而且是巴黎高等商学院（HEC Paris）的转学生。她的工作经验包括为摩根斯坦利从事结构化产品估值，同时承担波士顿咨询集团的管理顾问一职。她是一名活跃的外汇交易者，并且为Gazdaság电视台拍摄了一个月的投资思想的直播，在节目里她经常用到技术分析，这一主题自她15岁起就开始感兴趣。她曾经是维尔纽斯大学多门金融相关科目的助教。
Tamás Vadász拥有布达佩斯考文纽斯大学的经济学硕士学位。毕业之后，他从事金融服务业的顾问工作。目前，他正在进行金融学博士学位的学习，他的主要研究兴趣包括金融经济学和银行业的风险管理。他在考文纽斯大学教的课程包括金融计量学、投资学和公司财务。
Kata Váradi自2013年起是布达佩斯考文纽斯大学金融系的助理教授。作为金融学学生，Kata毕业于2009年，并在2012年其毕业论文关于匈牙利股票市场的市场流动性风险分析通过答辩，获得了布达佩斯考文纽斯大学的博士学位。她的研究领域是市场流动性、固定收益证券以及医疗保健系统的网络。除了做研究，她也积极从事教学。她主要讲授公司财务、投资学、估值以及跨国金融管理。
Ágnes Vidovics-Dancs是一位博士学位候选人，并且是布达佩斯考文纽斯大学的助理教授。此前她的工作是匈牙利政府债务管理局的初级风险管理师。她的主要研究领域是政府债务管理（一般）以及主权危机和违约（特别的）。她持有CEFA和CIIA证书。
Matthew Gilbert是一名量化分析师，在加拿大多伦多的加拿大养老基金投资公司（CPPIB）的全球宏观组工作。他拥有多伦多大学的量化金融硕士学位，以及皇后大学的应用数学和机械工程本科学位。
Hari Shanker Gupta博士是一位算法交易系统开发领域的高级量化研究分析师。在此之前，他是位于印度班加罗尔的印度科学研究院（Indian Institute of Science，IISc）的博士后。他在IISc获得了应用数学和科学计算的博士学位。他在印度瓦拉纳西的巴纳拉斯印度大学（Banaras Hindu University，BHU）完成了硕士学位。在瓦拉纳西，他因为杰出的表现，获得过4次金奖。
Hari已经在数学和科学计算的著名期刊上发表过5篇研究论文。他拥有在数学、统计和计算领域的工作经验。经验主题包括数值方法、偏微分方程、数理金融、随机分析、数据分析、时间序列分析、有限差分以及有限元方法。他擅长数学软件Matlab、统计编程语言R、Python和编程语言C。
他曾经为Packt出版的《量化金融R语言初级教程》（Introduction to R for Quantitative Finance）和《Python数据分析学习指南》（Learning Python Data Analysis）两本书审稿。
Ratan Mahanta拥有计算金融的硕士学位。目前他就职于GPSK投资集团，是一名高级量化分析师。对于卖方和风险顾问公司的量化交易和开发等方面，他拥有3年半的相关工作经验。他在Github开源平台为量化交易领域编写算法代码。他是自我激励、求知欲旺盛并且努力工作的人，热爱解决市场、技术、研究和设计交叉领域的困难问题。目前，他正在开发高频交易策略和量化交易策略。
量化交易：外汇、股权、期货与期权以及关于衍生品的工程。
算法：偏微分方程、随机微分方程、有限差分方法、蒙特卡洛以及机器学习。
代码：R编程、RStudio的Shiny、C++、Matlab、HPC以及科学计算。
数据分析：大数据分析（EOD到TBT）、Bloomberg、Quandl以及Quantopian。
策略：波动率套利、香草和奇异期权建模、趋势跟踪、均值回复、协整、蒙特卡洛模拟、在线价值、压力测试、高夏普比的买方策略、信用风险建模以及信用评分。
他也审阅过Packt出版社的《精通R语言科学计算》（Mastering Scientific Computing with R），目前，他正在审阅的图书是Packt出版社的《R语言机器学习手册》（Maching Learning with R Cookbook）。
在本章中，我们探讨一些时间序列分析的高级方法以及如何通过R来实现。作为一门学科，时间序列分析已有数百部著作，内容非常广泛（我们会在本章末的阅读列表中，列出在理论与R编程两方面最重要的参考目录）。我们责无旁贷地精心界定了本章的范围，专注于实证金融与量化交易中必不可少的重要主题。但是，在起始阶段我们必须强调，本章仅仅为时间序列分析的进一步研究奠定了基础。
我们之前曾经出版过一本图书——《量化金融R语言初级教程》（Introduction to R for Quantitative Finance）。那本书探讨了时间序列分析的一些基本主题，如线性单变量时间序列建模、自回归单整移动平均（Autoregressive Integrated Moving Average，ARIMA）建模、广义自回归条件异方差（Generalized Autoregressive Conditional Heteroskedasticity，GARCH）波动建模。如果你未曾用R做过时间序列分析，可以考虑同时阅读那本书的第1章“时间序列分析”。
本书对所有这些主题探讨得都更加深入。你将会熟悉一些重要概念，如协整、向量自回归模型、脉冲响应函数、非对称GARCH波动率建模（包括指数GARCH模型和门限GARCH模型）、信息冲击曲线。我们首先介绍相关理论，然后讲解多变量时间序列模型实际建模的知识，同时介绍几个有用的R包及其功能。此外，我们通过简单明了的例子，一步一步地指导读者使用R编程语言进行实证分析。
金融资产价格的运动、技术分析和量化交易的基本问题常常被纳入单变量框架下进行建模。我们能否预测证券价格未来是上升还是下降？这只特定的证券处于向上还是向下的趋势中？我们该买还是该卖？这些问题都需要慎重考虑。此外，投资者常常面对着更复杂的局面，不能仅仅把市场看成不相关的工具与决策问题组成的集合。
如果单独观察这些工具，可以发现正如市场有效假说所示，它们既非自相关又非均值可预测。但是，工具之间的相关性又显而易见。这个特性很可能为交易行为所利用，或者出于投机目的，或者出于对冲目的。以上探讨证实了多变量时间序列技术在数量金融中的应用。在本章中，我们会讨论两个在金融中应用广泛的著名计量经济学概念——协整和向量自回归模型。
现在，我们考虑一个时间序列向量
，它由元素
构成，每个元素独立地表示一个时间序列，如不同金融产品价格的变动。我们从协整数据序列的正式定义开始。
如果
的时间序列中的每个序列都是d阶单整（特别在大多数的应用中，这些序列是1阶单整，这意味着序列是非平稳的单位根过程或者随机游走过程）。同时，存在一个线性组合
是
阶单整（特别是当阶数为0时，序列是一个平稳过程）。那么，我们称时间序列向量为协整的。
直观地来看，这个定义意味着经济中存在某些潜藏力量。这种力量使n个时间序列看似是相互独立的随机游走过程，但在长期中趋于一致。协整时间序列的一个简单例子是下文这个取自Hamilton（1994）的配对向量。我们通过这个例子研究协整，同时熟悉R中一些基本模拟技术：
标准统计检验可以严格地展示中的单位根。R的tseries包或者urca包都可以实施单位根检验，本文使用后者。下文的R代码模拟了两个长度为1000的序列：
#生成两个长度为1000的时间序列
set.seed(20140623) #固定随机数种子
N <- 1000 #定义模拟的长度
x <- cumsum(rnorm(N)) #模拟一个正态随机游走
gamma <- 0.7 #设置初始的参数值
y <- gamma * x + rnorm(N) #模拟协整序列
plot(x, type='l') #画出两个序列
lines(y,col="red")
小提示
下载示例代码
如果你购买了Packt出版社的任何图书，都可以用自己的账号在网站http://www.packtpub.comhttp://www.packtpub.com下载示例的代码文件。如果你通过其他途径购买本书，那么可以访问网站http://www.packtpub.com/support，注册之后文件会直接通过邮件发送给你。
上列代码的输出结果如图1-1所示。
图1-1 模拟的协整序列
从图形上看，这两个序列看似都是随机游走过程，可以通过urca包实施ADF检验来检验平稳性。此外，R中还有很多其他检验方法。原假设的内容是过程中存在一个单位根（输出已经省略）。如果检验统计量小于临界值，我们就拒绝原假设：
#统计检验
install.packages('urca');library('urca')
#对模拟的单个时间序列做ADF test
summary(ur.df(x,type="none"))
summary(ur.df(y,type="none"))对于这两个模拟序列，检验统计量大于一般显著性水平（1%、5%和10%）的临界值。因此，我们不能拒绝原假设，并且我们认为这两个序列各自都是单位根过程。
现在，对这两个序列取下文中的线性组合并绘出残差序列：
z = y - gamma*x #取序列的线性组合
plot(z,type='l')上文代码的输出结果如图1-2所示。
图1-2 协整序列的线性组合
显而易见是个白噪声过程。但是拒绝单位根还是需要ADF检验的结果判定:
summary(ur.df(z,type="none"))在实际应用中，显然我们并不知道
的值。我们需要基于原始数据，通过一个序列对另一个序列进行线性回归来估计这个值。这个检验协整的方法称为Engle-Granger方法，实施这个方法需要下文两个步骤。
（1）对
进行线性回归（简单的最小二乘估计）。
（2）检验残差以确定单位根的存在。
小提示
我们在这里应该注意n个序列的情形，其中独立协整向量的可能个数落在0<r<n的范围。因此，当n>2时，协整关系可能不唯一。我们会在本章后面简要讨论r>1的情况。
简单线性回归可以通过lm函数拟合。残差可以从结果对象获得，如下文例子所示。ADF检验先按常规方式实施，再在所有显著性水平上确认是否拒绝原假设。此外，我们在本章后面讨论一些需要注意的事项：
#估计协整关系
coin <- lm(y ～ x -1) #不带截距项的回归
coin$resid #获取残差
summary(ur.df(coin$resid)) #残差的ADF检验现在，我们考虑如何将这个理论转化为成功的交易策略。此时，我们需要援引统计套利（statistical arbitrage）或者配对交易（pair trading）的概念。这个概念，在它最早期最简单的形式中，精准地探索这种协整关系。这些方法的主要目的是基于两个时间序列之间的价差来建立交易策略。如果序列是协整的，那么我们可以预期它们的平稳线性组合将会回复到0。因此，我们只要卖出高估证券买入低估证券，然后坐等价差回复，就能赚钱。
小提示
统计套利这个术语，通常用于多种复杂统计和计量经济学技术，目的是在统计意义上探索资产定价的相对偏误，而非在理论均衡模型中进行比较。
这种观点背后的经济原理是什么？潜在的经济力量决定了时间序列的线性组合形成协整关系。这种力量无法在统计模型中得到严格确认，但是有时会作为问题中变量之间的长期关系而存在。比如，我们可以预期相同行业中的相似公司会具有相似的增长，可以预期金融产品的即期价格和远期价格通过无套利原则联系在一起，还可以预期存在某种贸易往来的国家的外汇汇率会共同运动，或者可以预期短期利率和长期利率趋于接近。当交易员推测未来价格相关性的时候，变量偏离了统计上或者理论上的预期联动开启了各种各样的量化交易策略。
协整概念会在随后的章节中更深入讨论。为此，我们首先需要介绍向量自回归模型。
显然，向量自回归模型（VAR）可以看成自回归模型（AR）从单变量向多变量的扩展。它们在应用计量经济学中的普及可以追溯到Sims（1980）的开创性论文。VAR模型是最重要的多变量时间序列模型，在计量经济学和金融学中有广泛应用。R包中的vars为R用户提供了一个优秀的框架。如果需要这个包的详细评论，可以参考Pfaff（2013）。如果需要计量经济学理论，可以参考Hamilton（1994）、Lütkepohl（2007）、Tsay（2010）或者Martin et al.（2013）。在本书中，我们仅对该主题提供一个简洁又直观的概述。
在VAR模型中，我们从考虑一个长度为n的时间序列向量开始。VAR模型将每个变量的演进设定为所有其他变量滞后值的线性函数。换句话说，一个p阶的VAR模型如下：
这里，对所有的
，
是一个
的系数矩阵，而
是一个协方差矩阵正定的向量白噪声过程。向量白噪声这个术语假设不存在自相关，但允许各成分之间同期相关。换句话说，的协方差矩阵是非对角的。
矩阵记号清楚地解释了VAR模型的一个特殊的特征：所有变量仅仅取决于该变量本身与其他变量的滞后值，这意味着模型没有显式地刻画变量之间的同期相依性。该特征允许我们通过普通最小二乘方法逐个估计方程。这样的模型称为缩减型VAR模型，与之相对的是结构型VAR模型，我们随后介绍。
显然，不存在同期影响的假设存在过度简化的问题。而且，由此推出的脉冲—响应关系结果（某个特定变量发生冲击引起的过程变动）也许会产生误导性，也许没有用。这些问题引发了我们介绍结构型VAR模型，这一类模型显式地刻画了变量间的同期影响：
这里，
而且
，因此，结构形式可以通过简化形式乘上一个恰当的参数矩阵A得到，这个矩阵刻画了变量之间同期的结构关系。
小提示
和往常一样，我们的符号与vars包的技术文档保持一致，它与Lütkepohl（2007）的符号类似。
在缩减型模型中，同期相依性没有纳入建模。因此，这种相依性会出现在误差项的相关结构中，也就是的协方差矩阵，定义为
。在SVAR模型中，同期相依性显式地建模（通过左侧的矩阵A），并且扰动项定义为互不相关，因此协方差矩阵
是对角矩阵。在这里，扰动项常常指代结构冲击。
SVAR建模中既有趣又困难的地方是所谓的识别问题。SVAR模型是不可识别的，换句话说，如果没有额外约束，矩阵A的参数无法被估计。
小提示
我们应该如何理解模型是不可识别的？它的基本含义是指存在着不同的（无限多个）参数矩阵，却有相同的样本分布。因此，基于样本不能识别出参数的唯一值。
给定一个缩减式模型，总能导出一个恰当的参数矩阵，使残差正交。协方差矩阵
是半正定的，因此我们可以使用LDL分解（或者Cholesky分解也可以）。这种分解是指存在一个下三角矩阵L和一个对角阵D，使得
。通过选择
，结构模型的协方差矩阵变成
，由此给出
（注：原文此处有误）。现在，正如预期，我们得到
是一个对角阵的结论。注意，这种方法在本质上源于我们在方程上强加了一个主观的递归结构。irf()函数的缺省状态采用了这一方法。
在文献中，还有多种识别SVAR结构的方法，包括短期或者长期参数的约束，或者脉冲响应的符号约束［例如，见Fry-Pagan（2011）］。其中很多方法在R中仍然没有本地支持。这里，我们只介绍施加短期参数约束的一组标准技术，它们分别称为A模型、B模型和AB模型，每种技术在vars包里都有本地支持。
，对矩阵A上施加约束，使得
是一个对角协方差矩阵。为了保证模型是“恰好可识别”，我们需要
个额外约束。这让人想到施加了一个三角矩阵（但已不需要那种特定的结构）。
且
。这样基于约束模型残差有可能识别结构信息。
决定。脉冲响应分析常常是建立VAR模型的主要目标之一。脉冲响应函数本质上刻画了一个变量对系统中任何其他变量被冲击的反应（响应）。如果系统包含K个变量，那么可以确定
个脉冲响应函数。脉冲响应可以通过VAR过程的向量移动平均形式推导出数学方程，类似于单变量情形［细节参见Lutkepohl（2007）］。
作为一个说明性的案例，我们使用以下部分建立三元VAR模型。
我们的首要目的是通过额外变量预测股票市场指数，并识别脉冲响应。在这里，我们假设在给定的股票、股票市场整体和债券市场之间存在一种隐性长期关系。选择这个例子主要为了展示R编程环境中多种数据操作的可能性，以及通过一个很简单的案例说明详尽内容，而不是为了它的经济含义。
我们使用vars和quantmod包。如果你还没有这两个包，别忘了安装并加载它们：
install.packages('vars');library('vars')
install.packages('quantmod');library('quantmod')quantmod包提供了大量工具，可以直接从在线资源下载金融数据，在本书中我们会频繁地需要。我们使用getSymbols()函数：
getSymbols('MSFT', from='2004-01-02', to='2014-03-31')
getSymbols('SNP', from='2004-01-02', to='2014-03-31')
getSymbols('DTB3', src='FRED')默认状态下，雅虎财经（yahoofinance）用作股票和指数价格序列的数据源（参数设置为src='yahoo'，在本例中被省略了）。惯例是下载开盘价、最高价、最低价、收盘价、成交量及调整价格。下载的数据存储在一个xts数据类里，默认根据股票代码自动命名（MSFT和SNP）。通过调用泛型plot函数可以画出收盘价，但quantmod包的chartSeries函数提供了更好的图形示例。
使用以下缩写可以访问下载数据的各个部分：
Cl(MSFT) #收盘价
Op(MSFT) #开盘价
Hi(MSFT) #当日最高价
Lo(MSFT) #当日最低价
ClCl(MSFT) #收盘价对收盘价的日收益率
Ad(MSFT) #当日调整收盘价比如，使用这些缩写，日收盘价的收益率可以按以下方式绘制出来：
chartSeries(ClCl(MSFT)) #利用代码缩写的画图例子上述命令输出如图1-3所示。
图1-3 利用代码缩写MSFT绘制的股价图形
利率下载自联储经济数据（Federal Reserve Economic Data，FRED）数据源。接口的当前版本不允许对日期取子集。但是，下载的数据存储在一个xts数据类中，我们可以直接对需要的时期取子集。
DTB3.sub <- DTB3['2004-01-02/2014-03-31']下载的价格（假设是非平稳时间序列）需要先转化为平稳序列才能分析。换句话说，我们处理的是根据调整序列计算出的对数收益率：
MSFT.ret <- diff(log(Ad(MSFT)))
SNP.ret <- diff(log(Ad(SNP)))在转向VAR模型拟合之前，我们最后需要一个数据清洗的步骤才能继续。通过目测数据，我们可以看到在T-bill收益率序列中存在着缺失数据，而且数据集的长度各不相同（在某些日期，有利率报价，但没有股票价格）。为了解决这些数据质量问题，我们选择现在最易行的一种解决方案：合并数据集（通过省略那些没有全部3个数据的点），并删掉所有的NA数据。前者通过内连接参数执行（细节参见merge函数的帮助文档）。
dataDaily <- na.omit(merge(SNP.ret,MSFT.ret,DTB3.sub), join='inner')在这里需要注意，VAR建模通常对低频数据进行。有一种简单的方法，可以把数据转成为月度或者季度的频率，使用以下函数可以返回给定时期内的开盘价、最高价、最低价和收盘价：
SNP.M <- to.monthly(SNP.ret)$SNP.ret.Close
MSFT.M <- to.monthly(MSFT.ret)$MSFT.ret.Close
DTB3.M <- to.monthly(DTB3.sub)$DTB3.sub.Close简单的缩减VAR模型可以使用vars包的VAR()函数拟合数据。下列代码中的参数化允许在方程中有最大4阶的滞后，并且通过使用最优（最低）的赤池信息量值（Akaike Information Criterion）选择模型：
var1 <- VAR(dataDaily, lag.max=4, ic="AIC")对于更多已有模型的选择，可以考虑使用VARselect()，它提供了多种信息准则（输出省略）：
>VARselect(dataDaily,lag.max=4)结果对象是一个varest类的对象。通过summary()方法或者show()方法可以获得参数估计和其他更多统计结果（换句话说，只需输入变量）：
summary(var1)
var1还有其他方法值得一提。varest包的自定义绘图方法可以对所有变量——包括它的拟合值、残差、残差的自相关和偏自相关函数——分别生成图像。你需要敲击回车键获得新变量的图像。vars包还提供了大量自定义设置，请查阅这个包的文档。
plot(var1) #Diagram of fit and residuals for each variables
coef(var1) #concise summary of the estimated variables
residuals(var1) #list of residuals (of the corresponding ～lm)
fitted(var1) #list of fitted values
Phi(var1) #coefficient matrices of VMA representation通过简单地调用predict函数并加上想要的置信区间，可以得到使用估计的VAR模型进行的预测。
var.pred <- predict(var1, n.ahead=10, ci=0.95)脉冲响应由irf()函数首先生成数值结果，然后可以通过plot()方法绘出图形。接着我们可以再得到每个变量的不同图像，包括各个脉冲响应函数以及相应自助置信区间，如下列命令所示：
var.irf <- irf(var1)
plot(var.irf)现在，考虑使用之前讲过的A模型参数约束拟合结构VAR模型。识别SVAR模型需要
个约束，在我们的例子里，约束个数是3。
小提示
细节可以参见Lütkepohl（2007）。本来额外需要
个约束，但对角元素标准化为1，因而我们所需额外约束数如上。
SVAR模型的分析始于一个已估计的缩减型VAR模型（var1）。再用一个恰当的结构约束矩阵修正它。
为了简化，我们使用下列约束。
这些约束通过在矩阵A中设置相应的0，施加到SVAR模型中，矩阵A如下：
当在R中设定矩阵A为SVAR的估计参数，需要估计的参数的位因果设为NA值。这可以通过下列分配完成：
amat <- diag(3)
amat[2, 1] <- NA
amat[2, 3] <- NA
amat[3, 1] <- NA最后，我们可以拟合SVAR模型，并且绘出脉冲响应函数（输出省略）：
svar1 <- SVAR(var1, estmethod='direct', Amat = amat)
irf.svar1 <- irf(svar1)
plot(irf.svar1)最后，我们综合已学内容，并且讨论协整VAR和误差修正模型（Vector Error Correction Models，VECM）的概念。
我们从一个协整变量系统（例如，在交易情形下，这表示一组可能由相同基本面驱动的相似股票）开始。前面讨论过的标准VAR模型必须在变量平稳时才能估计。我们知道，移除单位根模型的传统方法是对序列取一阶差分。但是，在协整序列情形下，这会导致过度差分，并且损失了变量水平长期共同运动表达的信息。最终，我们的目的是建立一个平稳变量的模型，它也包含初始协整非平稳变量之间的长期联系，即建立一个协整VAR模型。向量误差修正模型（VECM）捕捉了这种思想，它包括变量差分的p-1阶VAR模型和一个从已知的（估计的）协整关系推导出的误差修正项。从直觉上，使用股票市场案例，VECM模型建立了股票收益率之间的短期关系，同时通过对价格长期共同变动的偏离加以修正。
一个两变量的VECM可以正式写成下面公式，并作为一个数值例子进行讨论。设
是一个由两个非平稳单位根序列
、
组成的向量，这两个序列协整，协整向量为
。因此，一个恰当的VECM模型可以如下面的公式所示：
这里，
且其中第一项通常称为误差修正项。
实践中有两种方法可以检验协整关系，并建立误差修正模型。在两变量情形下，Engle-Granger方法很有指导性，我们的数值案例基本沿袭了这种思想。在多变量情形下，其中可能的协整关系个数最大为(n-1)，这时需要采用Johansen过程进行协整检验。尽管后者的理论框架远超本书范围，我们仍然会简要介绍实用工具，并且为你更深的研究提供参考书。
为了讲授一些运用于VECM建模的基本R功能，我们使用一个标准案例。它包括3个月和6个月的国债二级市场利率，正如之前所讨论的，可以从FRED数据库下载。我们先任意选择某个时期，如从1984年～2014年，再考虑检验这段时间内的数据。ADF检验（Augmented Dickey Fuller tests）表明不能拒绝单位根原假设。
library('quantmod')
getSymbols('DTB3', src='FRED')
getSymbols('DTB6', src='FRED')
DTB3.sub = DTB3['1984-01-02/2014-03-31']
DTB6.sub = DTB6['1984-01-02/2014-03-31']
plot(DTB3.sub)
lines(DTB6.sub, col='red')我们通过运行一个简单的线性回归，可以一致地估计两个序列之间的协整关系。为了简化代码，我们定义变量
、
表示两个序列，y表示相应的向量序列。代码片段中的其他变量名惯例是明显的：
x1=as.numeric(na.omit(DTB3.sub))
x2=as.numeric(na.omit(DTB6.sub))
y = cbind(x1,x2)
cregr <- lm(x1 ～ x2)
r = cregr$residuals如果回归残差（变量r）恰好是变量的某个线性组合构成的平稳序列，那么这两个序列确实协整。你可以使用常见的ADF方法进行检验，但在这个情形下，传统的临界值不合适，应该使用修正值［例如，参见Phillips and Ouliaris（1990）］。
因此对于协整的存在性，更合适的方法是使用指定检验，如Phillips-Ouliaris检验，它在tseries和urca包中都可以运用。下面采用最基本的tseries包，演示如下：
install.packages('tseries');library('tseries');
po.coint <- po.test(y, demean = TRUE, lshort = TRUE)原假设表示两个序列不协整，因此p值小表示拒绝原假设并且协整关系存在。
当协整关系个数可能大于1时，适合使用Johansen过程。它的实施可以在urca包中找到：
yJoTest = ca.jo(y, type = c("trace"), ecdet = c("none"), K = 2)
######################
# Johansen-Procedure #
######################
Test type: trace statistic , with linear trend
Eigenvalues (lambda):
[1] 0.0160370678 0.0002322808
Values of teststatistic and critical values of test:
test 10pct 5pct 1pct
r <= 1 | 1.76 6.50 8.18 11.65
r = 0 | 124.00 15.66 17.95 23.52
Eigenvectors, normalised to first column:
(These are the cointegration relations)
DTB3.l2 DTB6.l2
DTB3.l2 1.000000 1.000000
DTB6.l2 -0.994407 -7.867356
Weights W:
(This is the loading matrix)
DTB3.l2 DTB6.l2
DTB3.d -0.037015853 3.079745e-05
DTB6.d -0.007297126 4.138248e-05对r=0（不存在协整关系）的检验统计量大于临界值，这表明拒绝原假设。然而，当
时，不能拒绝原假设。因此，我们认为存在协整关系。协整向量根据检验结果下方的标准化特征向量的第一列给出。
最后一步是得出这个系统的VECM表达式，换句话说，先计算协整方程，再求出差分变量的滞后和误差修正项，再对它们运行OLS回归。所选择的适合的函数里使用了我们之前创建的ca.jo对象类。下面代码里的参数r=1表示了协整的秩：
>yJoRegr = cajorls(dyTest, r=1)
>yJoRegr
$rlm
Call:
lm(formula = substitute(form1), data = data.mat)
Coefficients:
x1.d x2.d
ect1 -0.0370159 -0.0072971
constant -0.0041984 -0.0016892
x1.dl1 0.1277872 0.1538121
x2.dl1 0.0006551 -0.0390444
$beta
ect1
x1.l1 1.000000
x2.l1 -0.994407正如我们所预计的，误差修正项的系数是负的。从长期均衡水平上的短期的偏离会把变量推向零均衡偏离。
可以很容易地在二元情形下进行检查。Johansen过程会得到一个和采用了Engle-Grange过程从而一步步运用ECM相同的结果。这种观点在上传的R代码文件中有所示意。
金融时间序列的波动率会随着时间变化，这在实证金融中已经是广为熟悉和被接受的典型事实。但是，波动率的不可测性使得测量和预测它成为一项挑战性任务。通常，以下3种经验观察推动了波动率模型的演变。
波动性聚集：它指金融市场上的这样一种经验观察，平静期常常跟着平静期，而波动期常常跟着波动期。
资产收益率的非正态性：实证分析显示，相对于正态分布，资产收益率分布趋向于厚尾性。
杠杆效应：这会导致一种现象，波动率对正价格变动或负价格运动的反应往往不同。价格下降时波动率的增加幅度大于相似规模的价格上涨带来的波动率变动。
在下列代码中，我们演示了基于S&P资产价格的典型化事实。用我们已掌握的方法，从雅虎财经下载数据：
getSymbols("SNP", from="2004-01-01", to=Sys.Date())
chartSeries(Cl(SNP))我们的兴趣目标是日收益率序列，因此从收盘价计算对数收益率。其实quantmod包提供了一种更简单的方法，收益率可以直接计算：
ret <- dailyReturn(Cl(SNP), type='log')波动率分析始于观察自相关和偏自相关函数。我们希望对数收益率序列无关，但对数收益率的平方值或者绝对值都显示出了显著的自相关性。这意味着对数收益既不相关，也不独立。
注意下列代码中的par(mfrow=c(2,2))函数。通过它，可以重写R中默认的图形参数，把4张图形组织成方便的表格形式：
par(mfrow=c(2,2))
acf(ret, main="Return ACF");
pacf(ret, main="Return PACF");
acf(ret^2, main="Squared return ACF");
pacf(ret^2, main="Squared return PACF")
par(mfrow=c(1,1))上述代码输出如图1-4所示。
图1-4 收益率和平方收益率的ACF与PACF的比较
接下来，我们来看S&P的日度对数收益率的直方图和（或）经验分布，并与同均值同标准差的正态分布进行比较。我们使用函数density(ret)计算正态分布的非参数经验分布函数。再使用带有附加参数add=TRUE的函数curve()在刚刚画好的图形上绘出第二条线，如图1-5所示。
m=mean(ret);s=sd(ret);
par(mfrow=c(1,2))
hist(ret, nclass=40, freq=FALSE, main='Return histogram');curve(dnorm(x,
mean=m,sd=s), from = -0.3, to = 0.2, add=TRUE, col="red")
plot(density(ret), main='Return empirical distribution');curve(dnorm(x,
mean=m,sd=s), from = -0.3, to = 0.2, add=TRUE, col="red")
par(mfrow=c(1,1))图1-5 对数收益率的直方图和经验分布
很明显可以看出尖峰厚尾性。但我们还得在数值上确认（使用moments包）样本经验分布的峰度超过正态分布的峰度（等于3）。和其他的软件包不同，R报告名义峰度，而非超额峰度，结果如下：
> kurtosis(ret)
daily.returns
12.64 959放大图形的上尾和下尾可能有用，仅改变图形比例就可以，效果如图1-6所示。
# 放大尾部
plot(density(ret), main='Return EDF - upper tail', xlim = c(0.1, 0.2),
ylim=c(0,2));
curve(dnorm(x, mean=m,sd=s), from = -0.3, to = 0.2, add=TRUE, col="red")另一种有用的可视化练习是观看对数尺度上的密度（图1-7的左侧部分）或者Q-Q图（图1-7的右侧部分），它们是密度比较的常用工具。Q-Q图描绘了经验分位数对理论（正态）分布的图形。如果样本取自正态分布，它应该绘出一条直线。对这条直线的偏离表示存在厚尾性：
# 对数尺度上的密度图
plot(density(ret), xlim=c(-5*s,5*s),log='y', main='Density on log-scale')
curve(dnorm(x, mean=m,sd=s), from=-5*s, to=5*s, log="y", add=TRUE,
col="red")
# QQ-plot
qqnorm(ret);qqline(ret);图1-6 收益率经验分布的上尾
上述代码的输出如图1-7所示。
图1-7 对数尺度的收益率密度和Q-Q图
现在，我们可以把注意力转向波动率建模。
广义地说，金融计量经济学文献中有两类建模技术可以捕捉波动率的变化本质：GARCH族方法［Engle（1982）；Bollerslev（1986）］和随机波动模型（Stochastic Volatility，SV）。GARCH类模型和（纯粹的）SV类建模技术的主要区别在于，给定历史观测后，前者的条件方差是可以找到，而在SV模型中，考虑所有可得信息后，波动率依然无法预测。因此，SV模型的波动率本质上是隐性的，必须由测量方程滤出［比如，见Andersen–Benzoni（2011）］。换句话说，GARCH类模型意味着，可以通过历史观测值估计波动率，而在SV模型中，波动率自身有隐性的随机过程，因此推断潜在的波动率过程需要已实现收益率作为观测方程。
在本章中，基于两种原因我们介绍GARCH方法的基本建模技术。首先，它在应用工作中有更广泛的使用。其次，由于方法背景上的多样性，SV模型还不能被R包本地支持，而经验分析的使用需要大量的定制。
R提供了多个包用于GARCH建模，其中最著名的包是rugarch包、rmgarch包（处理多变量模型）以及fGarch包。但是，基本的tseries包也包含了一些GARCH功能。在本章中，我们会展示rugarch包的建模能力。本章中的标记方法遵从rugarch包中各自的输出和文档惯例。
GARCH(p,q)过程可以如下表示：
这里，
通常是条件均值方程（实际中常常是ARMA过程）的扰动项，并且
。换句话说，条件波动过程由它自己的滞后值
和滞后平方观测值（
的值）的线性组合决定。在实证研究中，GARCH（1,1）常常为数据提供了合适的拟合。我们将简单的GARCH（1,1）设定如下理解比较有益。模型中，条件方差通过长期方差
，最新一期预测方差
，以及新信息
的加权平均来设定［见Andersen et al.（2009）］。我们容易看出，GARCH模型如何捕捉到了波动的自回归特性（波动聚集性）和资产收益率分布的尖峰特性。但是，它的主要缺点在于它的对称性，因此无法捕捉分布的非对称特性与杠杆效应。
GARCH模型中出现波动聚集性高度符合直觉。
中的大的正（负）冲击会增加（减少）的值，进而会增加（减少）
的值，并导致更大（更小）值的
。冲击是持续的，因此这就是波动聚集。尖峰特性需要一定推导［见Tsay（2010）的例子］。
我们的实证案例是对苹果公司股票日收盘价收益率序列的分析，时间跨度从2006年1月1日开始，到2014年3月31日结束。在分析开始之前，作为一种有益的练习，我们建议你重复本章中的数据探索分析，以便确定苹果公司数据中的典型事实。
显然，第一步是安装包，如果它还未安装：
install.packages('rugarch');library('rugarch')通常，为了获取数据，我们通过quantmod包和getSymbols()函数，基于收盘价计算收益率序列：
#载入苹果股票数据并计算对数收益率
getSymbols("AAPL", from="2006-01-01", to="2014-03-31")
ret.aapl <- dailyReturn(Cl(AAPL), type='log')
chartSeries(ret.aapl)rugrach包的编程逻辑可以如下理解：无论你的目标是什么（拟合、滤波、预测还是模拟），你首先需要指定一个模型作为系统对象（变量），再依次将它插入各个函数。模型可以通过调用ugrachspec()函数设定。下文的代码设定了一个简单的GARCH（1,1）模型（sGARCH），均值方程中仅有一个常数
：
garch11.spec = ugarchspec(variance.model = list(model="sGARCH",
garchOrder=c(1,1)), mean.model = list(armaOrder=c(0,0)))一种很自然的处理方法是用日收益率时间序列数据拟合模型，即通过极大似然方法估计未知参数。
aapl.garch11.fit = ugarchfit(spec=garch11.spec, data=ret.aapl)这个函数在一系列输出中提供了
的参数估计：
> coef(aapl.garch11.fit)
mu omega alpha1 beta1
1.923328e-03 1.027753e-05 8.191681e-02 8.987108e-01我们可以通过生成对象（即仅仅通过键入这个变量名）的show()方法，获得估计和多种诊断检验。我们也可以通过键入适当的命令得到一大批的其他统计量、参数估计、标准误差和协方差矩阵估计。通过查阅ugarchfit对象类，可以得到完整的输出列表，下文代码展示了最重要的一部分：
coef(msft.garch11.fit) #估计的系数
vcov(msft.garch11.fit) #参数估计量的协方差矩阵
infocriteria(msft.garch11.fit) #常用信息量列表
newsimpact(msft.garch11.fit) #计算信息冲击曲线
signbias(msft.garch11.fit) #Engle - Ng符号偏差检验
fitted(msft.garch11.fit) #获得拟合的数据序列
residuals(msft.garch11.fit) #获得残差
uncvariance(msft.garch11.fit) #无条件的（长期）方差
uncmean(msft.garch11.fit) #无条件的（长期）均值标准的GARCH模型可以捕捉厚尾性和波动聚集性。但是，要想解释杠杆效应引起的非对称性，我们还需要更高级的模型。为了图示非对称性问题，我们接下来描述信息冲击曲线的概念。
信息冲击曲线，这一概念由Pagan和Schwert（1990）以及Engle和Ng（1991）提出，是一种图示波动对冲击反应的变化度量的有用工具。这个名字源于把新信息影响市场运动看成冲击的通常解释。它们画出了在不同规模冲击下条件波动的变化，可以简明地表达出波动的非对称影响。在如下的代码中，第一行对应着之前GARCH（1,1）模型定义计算的数值化的信息冲击，第二行创建了图形：
ni.garch11 <- newsimpact(aapl.garch11.fit)
plot(ni.garch11$zx, ni.garch11$zy, type="l", lwd=2, col="blue",
main="GARCH(1,1) - News Impact", ylab=ni.garch11$yexpr, xlab=ni.
garch11$xexpr)上述代码的输出截图如图1-8所示。
正如我们所料，无论对正冲击还是负冲击，波动的响应都不存在非对称性。现在，我们转向也能兼容非对称性的模型。
图1-8 由GARCH（1,1）模型计算的信息冲击曲线
Nelson（1991）提出了指数GARCH模型。这种方法直接对条件波动率的对数进行建模：
其中，E是期望算子。这个模型形式允许在变化的波动过程中存在乘法的动态。非对称性是由参数
刻画。负值表示过程对负冲击反应更大，正如在实际数据集的观察。
为了拟合EGARCH模型，模型设定中唯一需要改变的参数是设置EGARCH模型类型。通过运行fitting函数，可以估计出其他参数[见coef()]。
# 指定在均值方程中只带有一个常数的EGARCH(1,1) 模型
egarch11.spec = ugarchspec(variance.model = list(model="eGARCH",
garchOrder=c(1,1)), mean.model = list(armaOrder=c(0,0)))
aapl.egarch11.fit = ugarchfit(spec=egarch11.spec, data=ret.aapl)
> coef(aapl.egarch11.fit)
mu omega alpha1 beta1 gamma1
0.001446685 -0.291271433 -0.092855672 0.961968640 0.176796061如图1-9所示，信息冲击曲线反映了条件波动对于冲击反应的强烈非对称性，并且证实了非对称性模型的必要性。
ni.egarch11 <- newsimpact(aapl.egarch11.fit)
plot(ni.egarch11$zx, ni.egarch11$zy, type="l", lwd=2, col="blue",
main="EGARCH(1,1) - News Impact",
ylab=ni.egarch11$yexpr, xlab=ni.egarch11$xexpr)图1-9 由EGARCH（1,1）模型计算的信息冲击曲线
另一个著名的例子是TGARCH模型，解释更容易。TGARCH设定了模型参数在一个确定的门限之上和之下是不同的。TGARCH也是一种更一般的类（非对称的幂ARCH类）的子模型，因为它在应用金融计量经济学文献中的广泛深入的运用，我们单独地讨论它。
TGARCH模型的方程式如下：
其中
模型解释很直接。ARCH系数依赖于历史误差项的符号。如果
为正，负的误差项会对条件波动率有更高的影响，正如我们之前在杠杆效应所见。
在R的rugarch包中，门限GARCH模型在一类更一般的GARCH模型框架中应用，称为GARCH模型族［Ghalanos（2014）］。
# 指定在均值方程中只带有一个参数的TGARCH(1,1) 模型
tgarch11.spec = ugarchspec(variance.model = list(model="fGARCH",
submodel="TGARCH", garchOrder=c(1,1)),
mean.model = list(armaOrder=c(0,0)))
aapl.tgarch11.fit = ugarchfit(spec=tgarch11.spec, data=ret.aapl)
> coef(aapl.egarch11.fit)f
mu omega alpha1 beta1 gamma1
0.001446685 -0.291271433 -0.092855672 0.961968640 0.176796061由于特定的函数形式，门限GARCH的信息冲击曲线在表示不同响应时变化更小。运行下列命令，会看到在零点处有一个扭结，如图1-10所示。
ni.tgarch11 <- newsimpact(aapl.tgarch11.fit)
plot(ni.tgarch11$zx, ni.tgarch11$zy, type="l", lwd=2, col="blue",
main="TGARCH(1,1) - News Impact",
ylab=ni.tgarch11$yexpr, xlab=ni.tgarch11$xexpr)图1-10 由TGARCH（1,1）模型计算的信息冲击曲线
rugarch包允许从指定的模型以一种简单的方式模拟。当然，为了模拟，我们还需指定ugarchspec()内的模型参数，这可以通过fixed.pars参数完成。指定GARCH模型和给定条件均值之后，仅用ugarchpath()函数，就可以模拟出一个相应的时间序列：
garch11.spec = ugarchspec(variance.model = list(garchOrder=c(1,1)),
mean.model = list(armaOrder=c(0,0)),
fixed.pars=list(mu = 0, omega=0.1, alpha1=0.1,
beta1 = 0.7))
garch11.sim = ugarchpath(garch11.spec, n.sim=1000)一旦我们估计好严格拟合的模型，再来预测条件波动率只剩一步之遥：
aapl.garch11.fit = ugarchfit(spec=garch11.spec, data=ret.aapl, out.
sample=20)
aapl.garch11.fcst = ugarchforecast(aapl.garch11.fit, n.ahead=10,
n.roll=10)预测序列的绘图方法给用户提供了一个可选菜单，可以画出预测的时间序列或者预测的条件波动率，如图1-11所示。
plot(aapl.garch11.fcst, which='all')图1-11 预测的时间序列和条件波动率
在本章中，我们回顾了时间序列分析的某些重要内容，如协整、向量自回归和GARCH类条件方差模型。同时，我们介绍了一些R的有用知识和技巧，用来开始量化和实证金融的建模。我们希望你能从中获益，但是必须再次提到，不论是从时间序列和计量经济学理论的角度看，还是从R编程的角度看，本章内容皆不完善。互联网上有关于R编程语言的丰富文档资料，并且R的用户社区聚集着数以千计的高精尖人才。我们鼓励你能超越书本，成为一名自学者，遇到困难不要退缩。网络上肯定可以找到你处理困难需要的答案。要多多使用R包的文档和帮助文件，常常研究R的官网http://cran.r-project.org/。后续的章节会提供大量其他的示例，以帮助你在浩瀚的R功能、包和函数中找到自己的路。
金融资产的估值计算，大多数情况下基于现金流折现方法。即计算预期未来现金流的折现值，得到金融资产的现值。因此，为了能对资产估值，我们需要知道反映货币时间价值和给定资产风险的适当收益率。目前已有两种决定预期收益率的主流方法：资本资产定价模型（capital asset pricing model，CAPM）和套利定价理论（arbitrage pricing theory，APT）。CAPM是一个均衡模型，而APT建立在无套利原则之上。因而这两种方法的起点和内在逻辑都极为不同。但是，如果我们选择了合适的市场因素，通过两种方法所得到的最终定价公式却有可能极为相似。对于了解CAPM和APT的比较，可以参考Bodie-Kane-Marcus（2008）。当使用真实世界的数据来检验这些理论模型时，我们运行线性回归。由于我们在Daroczi et al.（2013）中已经详细讨论了CAPM，本章主要讨论APT。
本章分为两个部分。在前半部分中，我们先介绍一般的APT理论，然后再介绍一个特殊的三因素模型，这个三因素模型出自Fama和Frech的开创性论文。在后半部分中，我们说明了如何用R选择数据以及如何从实际市场数据中估计定价系数，最后我们使用更新的样本数据重新检验了著名的Fama-French模型。
APT基于这样的假设，市场中的资产收益率取决于宏观经济因素和公司特定因素。并且，资产收益率通过以下线性因素模型生成：
（方程1）
这里，
是资产i的预期收益率，
表示第j个因素的非预期变动，而
表示第i个证券对该因素的敏感性，同时
表示非预期的公司特定事件引起的收益。因此，
表示随机系统影响，表示非系统（即个体的）影响，非系统影响表示总影响中无法被系统因素捕捉的那部分。作为非预期的量，和都具有无条件的零均值。在这个模型中，包括系统特定风险在内的因素之间相互独立。因此，资产收益率分别来源于两部分：系统性风险因素影响了市场中的所有资产，非系统性风险仅仅影响特定公司。非系统性风险可以通过增加组合中资产种类来分散化。相反，来自经济整体的风险可以影响整个股票市场，系统性风险无法分散化（Brealey-Myers，2007）。
这个模型有这样一个推论，资产的已实现收益率是多个随机因素的线性组合（Wilmott，2007）。
APT的其他重要假设如下。
为1，而对所有其他因素的为0。从以上假设出发，可以推出任何组合的风险溢价等于因素组合的风险溢价加权和（Medvegyev-Szaz，2010）。下文的定价公式可以导出为两因素模型：
（方程2）
这里，
表示第i个资产的收益率，
表示无风险收益率，
表示第i个股票的风险溢价对第一个系统因素的敏感性，并且
表示这个因素的风险溢价。同样，
表示第i个股票的风险溢价对第二个因素的超额收益
的敏感性。
当我们需要实施APT时，使用如下形式的线性回归方程式：
（方程3）
这里，
表示一个常数，而i表示资产的非系统性的、公司特定的风险。所有其他变量含义如前所述。
如果模型中仅仅有一个因素，并且这个因素就是市场组合收益率，那么，CAPM模型和APT模型的定价方程式相同：
（方程4）
这种情形下，用真实市场数据检验的方程如下：
（方程5）
这里，
表示通过一个市场指数（如S&P500，即标准普尔500指数）代表的市场组合收益率。因此我们称方程5为指数模型。
APT的实现分4步进行：识别因素，估计因素系数，估计因素溢价，采用APT进行定价（Bodie et al.2008）。
（1）识别因素：因为APT本身不包含关于因素的任何内容，所以因素需要通过实证分析来识别。这些因素通常考虑宏观经济因素，如股票市场收益率、通货膨胀率、商业周期等。使用宏观经济因素的一大问题是各个因素相互不独立。因此常常需要使用因子分析来识别因素。但是，通过因子分析识别出的因素在经济学上不容易有好的解释。
（2）估计因素系数：为了估计多变量线性回归模型的系数，我们使用方程3的一个一般形式。
（3）估计因素溢价：因素溢价基于历史数据来估计，对因素组合溢价的历史时间序列数据取均值。
（4）给出APT定价方程：通过代入适合的变量，用方程2来计算任何资产的预期收益率。
Fama和French在1996年提出一个多因素模型。他们使用公司指标因素替代宏观因素，因为他们发现这些因素能够更好地描述资产的系统风险。Fama和French(1996)向市场组合收益率中增加了公司规模和净值市值比作为收益率生成因素，扩展了指数模型。
公司规模因素定义为小公司与大公司的收益率之差（
）。变量名是SMB，源于“small minus big”的首字母缩写。净值市值比因素定义为高净值市值比减去低净值市值比的公司收益率之差（
）。变量名是HML，源于“high minus low”的首字母缩写。
模型如下：
（方程6）
这里，是一个常数，表示异常收益率。是无风险收益率。
是第i个资产对净值市价比因素的敏感性。
是第i个股票的风险溢价对市场指数因素的敏感系数。
是市场指数因素的风险溢价。是资产的非系统性、公司特定风险，均值为零。
在接下来的部分中，我们将会学习在R的帮助下如何实现先前讲过的模型。
在第4章大数据—高级分析中，我们将会详细讨论获取开源数据及对其进行高效处理的各种知识与方法。在这里，我们仅仅讲述股票价格时间序列和其他相关信息如何获取和如何用于因素模型估计。
我们使用quantmod包来收集数据库。
以下是在R中实现的代码：
library(quantmod)
stocks <- stockSymbols()然后，我们需要等待几秒钟获取数据，接着可以看到输出：
Fetching AMEX symbols...
Fetching NASDAQ symbols...
Fetching NYSE symbols...现在，我们得到一个R数据框对象，这个对象涵盖了在各个交易所（如AMEX，NASDAQ，或者NYSE）进行交易的大约6500只股票。为了查看数据集涵盖的变量，我们用str命令：
str(stocks)
'data.frame': 6551 obs. of 8 variables:
$ Symbol : chr "AA-P" "AAMC" "AAU" "ACU" ...
$ Name : chr "Alcoa Inc." "Altisource Asset Management Corp"...
$ LastSale : num 87 1089.9 1.45 16.58 16.26 ...
$ MarketCap: num 0.00 2.44e+09 9.35e+07 5.33e+07 2.51e+07 ...
$ IPOyear : int NA NA NA 1988 NA NA NA NA NA NA ...
$ Sector : chr "Capital Goods" "Finance" "Basic Industries"...
$ Industry : chr "Metal Fabrications" "Real Estate"...
$ Exchange : chr "AMEX" "AMEX" "AMEX" "AMEX" ...我们可以删掉确实不需要的变量，也可以增加来自不同数据库的市场资本化信息和公司账面价值作为新变量，这些变量我们需要用来估计Fama-French模型：
stocks[1:5, c(1, 3:4, ncol(stocks))]
Symbol LastSale MarketCap BookValuePerShare
1 AA-P 87.30 0 0.03
2 AAMC 985.00 2207480545 -11.41
3 AAU 1.29 83209284 0.68
4 ACU 16.50 53003808 10.95
5 ACY 16.40 25309415 30.13我们也会需要无风险收益率的时间序列。我们通过计算月度LIBOR市场上的美元利率确定这个序列：
library(Quandl)
Warning message:
package 'Quandl' was built under R version 3.1.0
LIBOR <- Quandl('FED/RILSPDEPM01_N_B',
start_date = '2010-06-01', end_date = '2014-06-01')
Warning message:
In Quandl("FED/RILSPDEPM01_N_B", start_date = "2010-06-01", end_date
= "2014-06-01") : It would appear you aren't using an authentication
token. Please visit http://www.quandl.com/help/r or your usage may be
limited.因为数据依然分配给LIBOR变量，我们可以忽略这个警告信息。
Quandl包，tseries包和其他收集数据的包会在第4章（大数据—高级分析）中详细讨论。
这个方法也可以用来获取股票价格，并且标普500指数（S&P500指数）可以用来表示市场组合。
我们有一张表，记录了接近5000只股票价格在2010年6月1日到2014年6月1日之间的时间序列。其中第一列和最后几列如下所示：
d <- read.table("data.csv", header = TRUE, sep = ";")
d[1:7, c(1:5, (ncol(d) - 6):ncol(d))]
Date SP500 AAU ACU ACY ZMH ZNH ZOES ZQK ZTS ZX
1 2010.06.01 1070.71 0.96 11.30 20.64 54.17 21.55 NA 4.45 NA NA
2 2010.06.02 1098.38 0.95 11.70 20.85 55.10 21.79 NA 4.65 NA NA
3 2010.06.03 1102.83 0.97 11.86 20.90 55.23 21.63 NA 4.63 NA NA
4 2010.06.04 1064.88 0.93 11.65 18.95 53.18 20.88 NA 4.73 NA NA
5 2010.06.07 1050.47 0.97 11.45 19.03 52.66 20.24 NA 4.18 NA NA
6 2010.06.08 1062.00 0.98 11.35 18.25 52.99 20.96 NA 3.96 NA NA
7 2010.06.09 1055.69 0.98 11.90 18.35 53.22 20.45 NA 4.02 NA NA如果我们的数据已经储存在硬盘中，那么可以使用read.table函数方便的读取。在第4章（大数据—高级分析）中，我们会讨论如何从互联网收集数据。
现在，我们得到了所需数据：市场组合（标普500）、股票价格、无风险利率（月度LIBOR）。
为了清洗数据库，我们删除了有缺失值、0价格或者负价格的变量。最简单的实现如下：
d <- d[, colSums(is.na(d)) == 0]
d <- d[, c(T, colMins(d[, 2:ncol(d)]) > 0)]为了使用colMins函数，需要应用matrixStats包。现在，我们可以开始处理数据。
由于识别影响证券收益率的宏观变量很困难（Medvegyev-Szaz，2010，pp.42），我们在实践中使用因子分析相当不易。通常，我们通过主成分分析来寻找驱动收益率变动的潜因子。
在最初下载的6500只股票中，我们可以使用4500只股票数据。其他部分由于缺失值或者0价格的缘故删除了。现在，由于这里用不到日期，还有我们将标普500本身视为一个独立的因子，在主成分分析（principal component analysis，PCA）中不考虑，我们又删除了前两列数据。然后，计算对数收益率：
p <- d[, 3:ncol(d)]
r <- log(p[2:nrow(p), ] / p[1:(nrow(p) - 1), ])也有其他方法可以计算给定资产的对数收益率，就是PerformanceAnalytics库中的return.calculate(data,method="log")函数。
因为我们的股票数目过于庞大，为了实施PCA，要么需要数据时间长度至少25年，要么需要减少股票数量。使因素模型在几十年间保持稳定，这是不可能的。因此，为了达到说明的目的，我们随机选择了百分之十的股票，并对这个样本计算PCA模型：
r <- r[, runif(nrow(r)) < 0.1]runif(nrow(r)) < 0.1给出了一个4013维的0—1向量，选出了原表中几乎10%的列（在这个例子中，个数为393）。我们也可以使用以下样本函数，在网站http://stat.ethz.ch/R-manual/ R-devel/library/base/html/sample.html上你可以找到更多细节：
pca <- princomp(r)结果，我们得到一个princomp类对象。这个对象有8个属性，其中最重要的属性是加载矩阵和sdev属性（包括组成部分的标准差）。第一主成分是数据集方差最大的向量。
我们确认一下主成分的标准差：
plot(pca$sdev)结果如图2-1所示。
图2-1 主成分的标准差
我们可以看出，前5个成分是独立的。因此，我们应该选择5个因子。但是，其他因子的标准差同样显著。所以，不能通过几个因子解释整个市场。
我们可以通过调用factanal函数确认结果。这个函数估计了五因子模型：
factanal(r, 5)我们可以发现，实施这个计算花费了很多时间。因子分析与PCA相关，但从数学角度看稍复杂一些。结果，我们得到一个factanal类对象。它有很多属性。但是，此时我们仅对以下输出部分有兴趣：
Factor1 Factor2 Factor3 Factor4 Factor5
SS loadings 56.474 23.631 15.440 12.092 6.257
Proportion Var 0.144 0.060 0.039 0.031 0.016
Cumulative Var 0.144 0.204 0.243 0.274 0.290
Test of the hypothesis that 5 factors are sufficient.
The chi square statistic is 91756.72 on 75073 degrees of freedom.The
p-value is 0结果显示，五因子模型适合数据。但是，可解释的方差仅仅接近30%。这意味着模型应该考虑扩展进其他因子。
我们有一个包含4015只股票5年价格的数据框架，和包含LIBOR时间序列的LIBOR数据框架。首先，我们需要计算收益率，再与LIBOR利率合并。
第一步，我们删掉数学计算不需要的数据。然后，对保留下来的每一列计算对数收益率：
d2 <- d[, 2:ncol(d)]
d2 <- log(tail(d1, -1)/head(d1, -1))在计算了对数收益率之后，我们把日期放回到收益率中。然后，在最后一步合并两个数据集：
d <- cbind(d[2:nrow(d), 1], d2)
d <- merge(LIBOR, d, by = 1)需要提醒读者的是，对数据框架实施merge函数，相当于SQL中的（内）联语句。
结果如下：
print(d[1:5, 1:5])]
Date LIBOR SP500 AAU ACU
2010.06.02 0.4 0.025514387 -0.01047130 0.034786116
2010.06.03 0.4 0.004043236 0.02083409 0.013582552
2010.06.04 0.4 -0.035017487 -0.04211149 -0.017865214
2010.06.07 0.4 -0.013624434 0.04211149 -0.017316450
2010.06.08 0.4 0.010916240 0.01025650 -0.008771986我们将LIBOR利率转换为日度收益率：
d$LIBOR <- d$LIBOR / 36000由于LIBOR利率的报价基于货币市场基差——以及（实际天数/360）的天数计算约定——而且时间序列包含使用百分比表示的利率，我们把LIBOR除以36000。现在，我们需要计算Fama-French模型的3个变量。正如在数据选择部分中所讲，我们有股票的数据框：
d[1:5, c(1,(ncol(d) - 3):ncol(d))]
Symbol LastSale MarketCap BookValuePerShare
1 AA-P 87.30 0 0.03
2 AAMC 985.00 2207480545 -11.41
3 AAU 1.29 83209284 0.68
4 ACU 16.50 53003808 10.95
5 ACY 16.40 25309415 30.13我们得删掉那些没有价格数据的股票：
> stocks = stocks[stocks$Symbol %in% colnames(d),]我们将市场上限作为一个变量。我们仍需对每只股票计算账面市值比：
stocks$BookToMarketRatio <-
stocks$BookValuePerShare / stocks$LastSale
str(stocks)
'data.frame': 3982 obs. of 5 variables:
$ Symbol : Factor w/ 6551 levels "A","AA","AA-P",..: 14
72...
$ LastSale : num 1.29 16.5 16.4 2.32 4.05 ...
$ MarketCap : num 8.32e+07 5.30e+07 2.53e+07 1.16e+08...
$ BookValuePerShare: num 0.68 10.95 30.13 0.19 0.7 ...
$ BookToMarketRatio: num 0.5271 0.6636 1.8372 0.0819 0.1728 ...现在，我们需要计算SMB因素和HML因素。为了简化，我们将BIG公司定义为大于平均水平的公司。账面市值比实施同样的原则：
avg_size <- mean(stocks$MarketCap)
BIG <- as.character(stocks$Symbol[stocks$MarketCap > avg_size])
SMALL <- as.character(stocks[stocks$MarketCap < avg_size,1])这些数组包括了BIG公司和SMALL公司。现在，我们可以定义SMB因素：
d$SMB <- rowMeans(d[,colnames(d) %in% SMALL]) –
rowMeans(d[,colnames(d) %in% BIG])我们接着定义HML因素：
avg_btm <- mean(stocks$BookToMarketRatio)
HIGH <- as.character(
stocks[stocks$BookToMarketRatio > avg_btm, 1])
LOW <- as.character(
stocks[stocks$BookToMarketRatio < avg_btm, 1])
d$HML <- rowMeans(d[, colnames(d) %in% HIGH]) –
rowMeans(d[, colnames(d) %in% LOW])第三个因素这样计算：
d$Market <- d$SP500 - d$LIBOR定义完3个因素，我们在花旗集团（Citi）股票和伊克塞利克斯（EXEL）股票上试一试：
d$C <- d$C - d$LIBOR
model <- glm( formula = "C ～ Market + SMB + HML" , data = d)GLM（general linear model，一般线性模型）函数作用如下：它将数据和公式作为参数读入。公式是一个形为响应～条件的字符串，其中响应是数据框中的一个变量名，条件指定了模型中的预测子，它包含在数据集中通过操作符“＋”分隔开的变量名。这个函数也可以用于Logistic回归，只是缺省状态设定为线性。
模型输出如下：
Call: glm(formula = "C～Market+SMB+HML", data = d)
Coefficients:
(Intercept) Market SMB HML
0.001476 1.879100 0.401547 -0.263599
Degrees of Freedom: 1001 Total (i.e. Null); 998 Residual
Null Deviance: 5.74
Residual Deviance: 5.364 AIC: -2387模型概要输出如下：
summary(model)
Call:
glm(formula = "C～Market+SMB+HML", data = d)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.09344 -0.01104 -0.00289 0.00604 2.26882
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.001476 0.002321 0.636 0.525
Market 1.879100 0.231595 8.114 1.43e-15 ***
SMB 0.401547 0.670443 0.599 0.549
HML -0.263599 0.480205 -0.549 0.583
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.005374535)
Null deviance: 5.7397 on 1001 degrees of freedom
Residual deviance: 5.3638 on 998 degrees of freedom
AIC: -2387
Number of Fisher Scoring iterations: 2结果显示，唯一显著的因素是市场溢价。这表明花旗集团的股票收益率倾向于与整个市场本身共同变动。
使用以下命令可以画出结果：
estimation <- model$coefficients[1]+
model$coefficients[2] * d$Market +
model$coefficients[3]*d$SMB +
model$coefficients[4]*d$HML
plot(estimation, d$C, xlab = "estimated risk-premium",
ylab = "observed riks premium",
main = "Fama-French model for Citigroup")
lines(c(-1, 1), c(-1, 1), col = "red")图2-2绘出了对花旗集团的Fama-French模型估计的风险溢价。
图2-2 花旗集团的Fama-French模型
看图2-2可以发现，收益率中存在一个异常值。我们将这个异常值设为0，看看不考虑它之后的结果：
outlier <- which.max(d$C)
d$C[outlier] <- 0如果我们运行相同的代码建立模型，并再次计算收益率的估计值和观测值，得到以下结果：
model_new <- glm( formula = "C ～ Market + SMB + HML" , data = d)
summary(model_new)
Call:
glm(formula = "C ～ Market + SMB + HML", data = d)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.091733 -0.007827 -0.000633 0.007972 0.075853
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0000864 0.0004498 -0.192 0.847703
Market 2.0726607 0.0526659 39.355 < 2e-16 ***
SMB 0.4275055 0.1252917 3.412 0.000671 ***
HML 1.7601956 0.2031631 8.664 < 2e-16 ***
---
Signif. codes: 0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘ ’1
(Dispersion parameter for gaussian family taken to be 0.0001955113)
Null deviance: 0.55073 on 1001 degrees of freedom
Residual deviance: 0.19512 on 998 degrees of freedom
AIC: -5707.4
Number of Fisher Scoring iterations: 2根据以上结果，所有3个因素均显著。
GLM函数不返回
。lm函数对线性回归同样可以得到精确值。我们可以从模型总结中读出r.squared = 0.6446。
结果显示，变量可以解释花旗集团风险溢价中超过64%的变动。我们绘出新结果：
estimation_new <- model_new$coefficients[1]+
model_new$coefficients[2] * d$Market +
model_new$coefficients[3]*d$SMB +
model_new$coefficients[4]*d$HML
dev.new()
plot(estimation_new, d$C, xlab = "estimated risk-premium",ylab =
"observed riks premium",main = "Fama-French model for Citigroup")
lines(c(-1, 1), c(-1, 1), col = "red")这个例子的输出如图2-3所示。
图2-3 花旗集团去掉异常值后的Fama-French模型
我们再检验另一只股票EXEL：
d$EXEL <- d$EXEL–d$LIBOR
model2 <- glm( formula = "EXEL～Market+SMB+HML" , data = d)
Call: glm(formula = "EXEL～Market+SMB+HML", data = d)
Coefficients:
(Intercept) Market SMB HML
-0.001048 2.038001 2.807804 -0.354592
Degrees of Freedom: 1001 Total (i.e. Null); 998 Residual
Null Deviance: 1.868
Residual Deviance: 1.364 AIC: -3759模型小结输出如下：
summary(model2)
Call:
glm(formula = "EXEL～Market+SMB+HML", data = d)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.47367 -0.01480 -0.00088 0.01500 0.25348
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.001773 0.001185 -1.495 0.13515
Market 1.843306 0.138801 13.280 < 2e-16 ***
SMB 2.939550 0.330207 8.902 < 2e-16 ***
HML -1.603046 0.535437 -2.994 0.00282 **
---
Signif. codes: 0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘ ’1
(Dispersion parameter for gaussian family taken to be 0.001357998)
Null deviance: 1.8681 on 1001 degrees of freedom
Residual deviance: 1.3553 on 998 degrees of freedom
AIC: -3765.4
Number of Fisher Scoring iterations: 2根据以上结果，所有3个因子均显著。
GLM函数并不包括。但lm函数对线性模型也可以得到精确的结果。我们从模型小结中得到r.squared = 0.2723。根据这个结果，我们认为变量可以解释EXEL风险溢价的超过27%的变动。
使用以下命令可以绘出图2-4：
estimation2 <- model2$coefficients[1] +
model2$coefficients[2] * d$Market +
model2$coefficients[3] * d$SMB + model2$coefficients[4] * d$HML
plot(estimation2, d$EXEL, xlab = "estimated risk-premium",
ylab = "observed riks premium",
main = "Fama-French model for EXEL")
lines(c(-1, 1), c(-1, 1), col = "red")图2-4 EXEL的Fama-French模型
本章中，我们看到如何建立并实现多因素模型。通过主成分分析，我们确认需要5个独立因素解释资产定价。但是，这些因素仅有30%的方差解释力，表现出模型解释的不充分性。通过实例，我们用实际市场数据重建了著名的Fama-French模型。在这个模型中，除了市场因素，我们也使用了另两个公司特定因素（SMB和HML）。我们使用内置函数进行主成分分析和因子分析，并讲述了如何使用一般线性模型进行回归分析。
我们发现，这3个因素是显著的。因此我们得出结论，对新近的样本，Fama-French因素具有解释力。我们非常鼓励你能够模仿经典理论，发展并检验新的多因素定价方程式，可能会比经典方程更加出色。
