原文：
zh.annas-archive.org/md5/f22c78458c840dd6848aa3910b34a7da译者：飞龙
前言
量化金融中的 R 入门 将展示如何使用统计计算语言 R 和 QF 来解决实际的量化金融问题。本书将涉及从时间序列分析到金融网络的各类话题。每一章将简要介绍相关理论，并通过 R 来解决特定问题。
本书内容
第一章，时间序列分析（Michael Puhle），解释了如何在 R 中处理时间序列数据。此外，你还将学习如何建模和预测房价，利用协整提高对冲比率，并建模波动性。
第二章，投资组合优化（Péter Csóka, Ferenc Illés, Gergely Daróczi），介绍了投资组合选择背后的理论思想，并展示了如何将这些知识应用于实际数据。
第三章，资产定价模型（Kata Váradi, Barbara Mária Dömötör, Gergely Daróczi），在前一章的基础上，介绍了资产回报与风险之间关系的模型。我们将讨论资本资产定价模型和套利定价理论。
第四章，固定收益证券（Márton Michaletzky, Gergely Daróczi），讲解了固定收益工具的基础知识。此外，你还将学习如何计算这种工具的风险，并构建能够抵御利率变化的投资组合。
第五章，利率期限结构估计（Tamás Makara, Gergely Daróczi），介绍了收益率曲线的概念，并展示了如何利用政府债券的价格来估算该曲线。
第六章，衍生品定价（Ágnes Vidovics-Dancs, Gergely Daróczi），解释了如何使用离散时间和连续时间模型进行衍生品定价。此外，你将学习如何计算衍生品的风险度量以及所谓的“希腊字母”。
第七章，信用风险管理（Dániel Havran, Gergely Daróczi），介绍了信用违约模型，并展示了如何使用 Copula 模型来建模相关的违约事件。
第八章，极值理论（Zsolt Tulassay），介绍了极值理论在保险和金融中的潜在应用。你将学习如何为火灾损失分布的尾部拟合模型。然后，我们将使用拟合的模型来计算风险价值（Value-at-Risk）和预期亏损（Expected Shortfall）。
第九章，金融网络（Edina Berlinger, Gergely Daróczi）解释了如何在 R 中表示、模拟、可视化和分析金融网络。我们将分析银行间借贷市场，并学习如何系统地识别重要的金融机构。
本书所需的内容
本书提供的所有代码示例应在安装好的 R 控制台中运行。你可以免费下载该软件，并在 r-project.org 上找到所有主要操作系统的安装说明。尽管我们不会涉及如何在集成开发环境中使用 R 等高级主题，但除了其他编辑器外，Emacs、Eclipse、vi 或 Notepad++ 都有很棒的插件，我们也强烈推荐尝试 RStudio，它是一个专为 R 设计的免费开源 IDE。
除了一个正常工作的 R 安装外，我们还将使用一些用户贡献的 R 包，这些包可以从综合 R 包库轻松安装。要安装一个包，可以在 R 控制台中使用 install.packages 命令，示例如下：
> install.packages('zoo')
安装之后，该包在使用前也应先加载到当前的 R 会话中：
> library(zoo)
你可以在 R 的主页上找到免费的入门文章和手册，但本书是面向初学者的，因此假设读者没有额外的 R 知识。
本书适合谁阅读
本书面向希望使用 R 解决定量金融问题的读者。假设读者对金融有所了解，但我们通常也会提供金融理论知识。假设读者不熟悉 R。希望开始学习 R 的读者可能会发现本书非常有用，因为我们不会全面介绍 R 语言，而是展示如何使用其中的一部分来解决特定问题。即使你已经在使用 R，你也会惊讶于它能应用于如此广泛的问题。
约定
本书中，你将会看到多种文本样式，用以区分不同类型的信息。以下是这些样式的示例，并解释其含义。
文本中的代码词汇、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟网址、用户输入和 Twitter 账号会如下所示：“我们将使用 forecast 包中的一些方法”
一段 R 代码（通常是函数体）如下所示：
logreturn <- function(x) {
log(tail(x, -1) / head(x, -1))
}
当我们希望引起你对代码块中特定部分的注意时，相关行或项会以粗体显示：
logreturn <- function(x) {
log(tail(x, -1) / head(x, -1))
}
任何命令行输入或输出都会按照以下方式编写：
> pi
[1] 3.141593
其中 “>” 表示 R 控制台正在等待命令的执行。多行表达式在第一行以相同的符号开始，但其余所有行的开头都会有 “+” 符号，以表示上一个 R 表达式仍未完成。
新术语 和 重要词汇 会以粗体显示。你在屏幕上看到的单词，比如在菜单或对话框中出现的单词，将以这种方式出现在文本中：“点击 下一步 按钮会将你带到下一个屏幕”。
注意
警告或重要提示会以框的形式呈现。
提示
提示和技巧将以这种方式呈现。
读者反馈
我们始终欢迎读者的反馈。告诉我们你对本书的看法——你喜欢的地方或可能不喜欢的地方。读者反馈对我们开发真正能让你受益的书籍至关重要。
要向我们发送一般反馈，只需发送电子邮件至 <feedback@packtpub.com>，并在邮件主题中提及书名。
如果你在某个领域拥有专业知识，并且有兴趣写作或为书籍做贡献，请参考我们的作者指南，访问 www.packtpub.com/authors。
客户支持
既然你已经是一本 Packt 书籍的骄傲拥有者，我们有很多资源帮助你充分利用这次购买。
下载示例代码
你可以从你的账户中下载你所购买的所有 Packt 书籍的示例代码文件，访问 www.packtpub.com。如果你是在其他地方购买了本书，可以访问 www.packtpub.com/support 注册并直接将文件通过电子邮件发送给你。
勘误
虽然我们已经非常小心地确保内容的准确性，但错误仍然会发生。如果你在我们的书籍中发现错误——可能是文本或代码中的错误——我们将非常感激你能向我们报告。这样，你不仅能帮助其他读者避免困扰，还能帮助我们改进后续版本的书籍。如果你发现任何勘误，请通过访问 www.packtpub.com/submit-errata 来报告，选择你的书籍，点击 勘误提交表单 链接，并输入你的勘误详情。勘误经过验证后，我们会接受并将其上传至我们的网站，或者将其添加到该书标题下的现有勘误列表中。你可以通过选择书名，在 www.packtpub.com/support 查看任何现有的勘误。
盗版
网络上的版权材料盗版问题在所有媒体中都普遍存在。在 Packt，我们非常重视保护我们的版权和许可。如果你在互联网上发现任何非法复制的我们作品，请立即提供相关的地址或网站名称，以便我们采取措施。
请通过 <copyright@packtpub.com> 联系我们，并附上涉嫌盗版材料的链接。
我们感谢你在保护我们的作者以及我们为你提供有价值内容的能力方面的帮助。
问题
如果您在书籍的任何方面遇到问题，可以通过<questions@packtpub.com>与我们联系，我们将尽力解决。
第一章：时间序列分析
时间序列分析关注的是随着时间推移收集的数据的分析。相邻的观察值通常是相互依赖的。因此，时间序列分析处理的是分析这种依赖关系的技术。
本章的目标是通过具体应用介绍一些常见的建模技术。我们将看到如何使用 R 来解决这些现实世界的例子。我们从一些关于如何在 R 中存储和处理时间序列数据的思考开始。接下来，我们将处理线性时间序列分析，并讨论如何利用它来建模和预测房价。在随后的章节中，我们利用协整的概念，通过考虑长期趋势，改进基本的最小方差对冲比率。最后，章节以如何使用波动率模型进行风险管理的部分作为结尾。
处理时间序列数据
原生的 R 类适用于存储时间序列数据，包括 vector、matrix、data.frame 和 ts 对象。但这些对象可以存储的数据类型有限，此外这些表示方法提供的功能也有限。幸运的是，存在一些专门处理时间序列数据更广泛表示的对象：zoo、xts 或 timeSeries 对象，这些都可以从同名包中获得。
并不一定每一个时间序列分析问题都需要创建时间序列对象，但更复杂的分析确实需要时间序列对象。你可以计算作为 R 向量表示的时间序列数据的均值或方差，但如果你想使用 decompose 进行季节性分解，你就需要将数据存储在时间序列对象中。
在以下示例中，我们假设你正在使用 zoo 对象，因为我们认为它是最广泛使用的包之一。在使用 zoo 对象之前，我们需要安装并加载 zoo 包（如果你已经安装了它，只需要加载即可），可以使用以下命令：
> install.packages("zoo")
> library("zoo")
为了熟悉可用的方法，我们从 Apple 股票的每日收盘价创建一个名为 aapl 的 zoo 对象，这些数据存储在 CSV 文件 aapl.csv 中。每一行包含一个日期和一个收盘价，两者之间用逗号分隔。第一行包含列标题（Date 和 Close）。日期的格式按照 ISO 8601 的推荐标准（YYYY-MM-DD）进行。收盘价已经调整了股票分拆、股息和相关变化。
提示
下载示例代码
你可以从自己账户中的www.packtpub.com下载所有已购买 Packt 书籍的示例代码文件。如果你在其他地方购买了本书，可以访问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="%")
我们可以将分析限制在时间序列的一个子集（一个窗口）中。苹果公司在 2013 年的最高股价可以通过以下命令行来查找：
> aapl_2013 <- window(aapl, start = '2013-01-01', end = '2013-12-31')
> aapl_2013[which.max(aapl_2013)]
2013-01-02
545.85
从风险管理的角度来看，回报分布的分位数是非常重要的。例如，我们可以通过一种简单的历史方法轻松确定 1 天 99%的风险价值。
> quantile(ret_simple, probs = 0.01)
1%
-7.042678
因此，任何一天的回报低于 7%的概率只有 1%。但是，如果这一天确实发生（大约每年会发生 2.5 次），那么你将至少损失 7%。
线性时间序列建模与预测
一类重要的线性时间序列模型是**自回归积分滑动平均（ARIMA）**模型，由Box 和 *Jenkins (1976)*提出。它假设当前的值仅能依赖于时间序列自身的过去值，或者依赖于某些误差项的过去值。
根据 Box 和 Jenkins，建立 ARIMA 模型包括三个阶段：
-
模型识别。
-
模型估计。
-
模型诊断检查。
模型识别步骤涉及使用图形方法或信息准则来确定一个初步模型的阶数（包含的过去值的数量以及过去误差项的数量）。在确定模型的阶数之后，模型的参数需要进行估计，通常使用最小二乘法或最大似然法。然后，必须仔细检查拟合的模型，以确保不存在模型不适合的情况。这是通过确保模型残差表现为白噪声来完成的；也就是说，残差中不再存在线性依赖。
英国房价建模与预测
除了zoo包，我们还将使用forecast包中的一些方法。如果你还没有安装它，你需要使用以下命令进行安装：
> install.packages("forecast")
之后，我们需要使用以下命令加载类：
> library("forecast")
首先，我们将每月的房价数据（来源：全国建筑协会）存储在一个zoo时间序列对象中。
> hp <- read.zoo("UKHP.csv", sep = ",",+ header = TRUE, format = "%Y-%m", FUN = as.yearmon)
FUN参数将给定的函数（as.yearmon，表示每月数据点）应用到日期列。为了确保我们确实存储了每月数据（每个周期有 12 个子周期），通过指定as.yearmon，我们查询数据序列的频率。
> frequency(hp)
[1] 12
结果意味着我们在一个周期（称为“年”）中有十二个子周期（称为“月”）。我们再次使用简单收益率进行分析。
> hp_ret <- diff(hp) / lag(hp, k = -1) * 100
模型识别与估计
我们使用forecast包提供的auto.arima函数来一次性识别最佳模型并估计系数。除了收益率序列（hp_ret），该函数还接受多个参数。通过指定stationary = TRUE，我们将搜索限制为平稳模型。类似地，seasonal = FALSE将搜索限制为非季节性模型。此外，我们选择 Akaike 信息准则作为相对质量的度量来用于模型选择。
> 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² estimated as 1.105: log likelihood=-390.97
AIC=789.94 AICc=790.1 BIC=804.28
根据赤池信息量准则（AIC），一个 AR(2)过程似乎最适合这些数据。为了进行视觉验证，我们可以使用命令pacf绘制偏自相关函数。该图表明直到滞后二，偏自相关值非零，因此，二阶 AR 过程似乎是合适的。这两个 AR 系数、截距（如果模型包含 AR 项的话，实际上是均值）以及各自的标准误差都已给出。在以下示例中，它们在 5%的显著性水平下均显著，因为相应的置信区间不包含零：
> confint(mod)
2.5 % 97.5 %
ar1 0.1174881 0.3422486
ar2 0.2364347 0.4617421
intercept 0.1368785 0.7321623
如果模型包含不显著的系数，我们可以使用arima函数重新估计模型，并使用fixed参数输入一个包含0和NA的元素向量。NA表示相应的系数应当被估计，而0表示该系数应设置为零。
模型诊断检查
快速验证模型的一种方法是使用以下命令绘制时间序列诊断图：
> tsdiag(mod)
前一个命令的输出如图所示：
我们的模型表现良好，因为标准化残差没有表现出波动性聚集，残差之间没有显著的自相关（根据自相关图 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")
输出结果如下图所示：
此外，我们还可以计算常见的准确性度量。
> accuracy(mod)
ME RMSE MAE MPE MAPE MASE
0.00120 1.0514 0.8059 -Inf Inf 0.792980241
该命令返回平均误差、均方根误差、平均绝对误差、平均百分比误差、平均绝对百分比误差和平均绝对标准化误差。
预测
要预测未来三个月（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
因此，我们预计未来三个月的平均房价会略有上涨，但标准误差较高，大约为 1%。为了绘制带标准误差的预测图，我们可以使用以下命令：
> plot(forecast(mod))
协整
协整的概念由*Granger（1981）提出，并由Engle 和 Granger（1987）*正式化，其核心思想是找到非平稳时间序列之间的线性组合，从而得到一个平稳的时间序列。因此，可以检测非平稳时间序列（例如价格）之间的稳定长期关系。
跨期对冲喷气燃料
航空公司是喷气燃料的天然买家。由于喷气燃料价格可能非常波动，大多数航空公司都会对其喷气燃料价格变化的部分风险进行对冲。在缺乏液态喷气燃料场外交易工具的情况下，航空公司使用相关的交易所交易期货合约（例如取暖油）进行对冲。在接下来的部分中，我们首先采用经典方法，仅考虑两种价格之间的短期波动，推导出最佳对冲比率；然后，我们通过考虑两者价格之间的长期稳定关系，改进经典的对冲比率。
我们首先加载必要的库。urca库提供了一些用于单位根检验和估计协整关系的有用方法。
> library("zoo")
> install.packages("urca")
> library("urca")
我们导入了喷气燃料和取暖油的月度价格数据（以每加仑美元计）。
> prices <- read.zoo("JetFuelHedging.csv", sep = ",",+ FUN = as.yearmon, format = "%Y-%m", header = TRUE)
仅考虑两种商品的短期行为（每月价格变化），可以通过拟合一个线性模型，解释喷气燃料价格变化由取暖油价格变化引起，从而推导出最小方差对冲。该回归的β系数即为最佳对冲比率。
> simple_mod <- lm(diff(prices$JetFuel) ~ diff(prices$HeatingOil)+0)
lm函数（用于线性模型）估计喷气燃料价格变化与取暖油价格变化的最佳拟合系数。+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")
我们使用恩格尔和格兰杰的两步估计技术。首先，使用扩展的迪基-富勒检验对两个时间序列进行单位根（非平稳性）检验。
> 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%的显著性水平下无法被拒绝，因为检验统计量-1.1335 并不比临界值-3.46 更负。对取暖油价格而言也是如此（检验统计量为-1.041）。
> ho_adf <- ur.df(prices$HeatingOil, type = "drift")
> summary(ho_adf)
我们现在可以继续估计静态均衡模型，并使用扩展的迪基-富勒检验对残差进行平稳性检验。请注意，由于所研究的序列是估计值，因此现在必须使用不同的临界值[例如，来自Engle 和 Yoo (1987)]。
> mod_static <- summary(lm(prices$JetFuel ~ prices$HeatingOil))
> error <- residuals(mod_static)
> error_cadf <- ur.df(error, type = "none")
> summary(error_cadf)
获得的检验统计量为-8.912，而在 1%水平下，样本大小为 200 时的临界值为-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）：价格之间的大幅偏离将被修正，价格将趋向于它们的长期稳定关系。
建立波动率模型
正如我们之前所见，ARIMA 模型用于对一个过程的条件期望进行建模，前提是已知其过去。对于这种过程，条件方差是常数。现实世界中的金融时间序列呈现出波动聚集现象；也就是说，相对平静的时期会被波动的爆发所打断。
本节我们将探讨 GARCH 时间序列模型，这些模型可以考虑到现实世界（金融）时间序列的这一特征，并将这些模型应用于 VaR 预测。
风险管理中的波动率预测
金融机构通过使用 VaR（在险价值）来衡量其活动的风险，通常在 99%的置信水平下计算，并以 10 个交易日的期限进行估算。这是预期仅在 1%的时间内超过的损失。
我们加载zoo库并导入从 1973 年 1 月到 2008 年 12 月的英特尔公司月度收益数据。
> library("zoo")
> intc <- read.zoo("intc.csv", header = TRUE,+ sep = ",", format = "%Y-%m", FUN = as.yearmon)
ARCH 效应检验
收益的图表表明，月度收益数据中可能存在 ARCH 效应。
> plot(intc, main = "Monthly returns of Intel Corporation",+ xlab = "Date", ylab = "Return in percent")
前面命令的输出如以下图所示：
我们可以使用统计假设检验来验证我们的直觉。两种常用的检验如下：
-
用于平方收益的自相关的 Ljung-Box 检验（作为波动率的代理）
-
拉格朗日乘数（LM）检验由*Engle（1982）*提出
首先，我们使用以下命令对平方收益的前 12 个滞后期进行 Ljung-Box 检验：
> Box.test(coredata(intc²), type = "Ljung-Box", lag = 12)
Box-Ljung test
data: coredata(intc²)
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 模型。
GARCH 模型规格
最常用的 GARCH 模型，也是通常适用于金融时间序列的模型，是 GARCH(1,1)模型。我们使用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)))
GARCH 模型估计
通过最大似然法对系数进行实际拟合的函数是ugarchfit，它以模型规格和回报数据作为输入。
> intc_garch11_fit <- ugarchfit(spec = intc_garch11_spec,+ data = intc)
关于更多参数，请参见ugarchfit的帮助文档。拟合模型的输出（使用命令intc_garch11_fit）提供了有用的信息，例如最优参数的值、对数似然函数值和信息准则。
风险模型回测
检查模型性能的有用测试是进行历史回测。在风险模型回测中，我们将估算的 VaR 与该期间的实际回报进行比较。如果回报比 VaR 更为负值，则说明 VaR 超出。在我们的案例中，VaR 超出应只发生在 1%的情况下（因为我们指定了 99%的置信水平）。
函数ugarchroll对指定的 GARCH 模型（这里的模型是intc_garch11_spec）执行历史回测。我们按照如下方式指定回测：
-
将使用的回报数据存储在
zoo对象intc中 -
回测的起始期（
n.start）应为序列开始后的 120 个月（即 1983 年 1 月） -
模型应每月重新估计（
refit.every = 1） -
我们使用
moving窗口进行估计 -
我们使用
hybrid求解器 -
我们希望计算 99%VaR 尾部水平的 VaR（
calculate.VaR = TRUE）（VaR.alpha = 0.01） -
我们希望保留估计的系数（
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: 5Actual %: 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 检验是对无条件覆盖率和超出值独立性的联合检验。在我们的案例中，尽管实际有五次超出，而预期为三次，我们仍不能拒绝零假设，即超出值是正确且独立的。
回测性能的图形也可以轻松生成。首先，使用从ugarchroll对象中提取的预测 VaR 创建一个zoo对象。
> intc_VaR <- zoo(intc_garch11_roll@forecast$VaR[, 1])
我们还会将zoo对象的index属性覆盖为该对象的rownames（年和月）。
> 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))
现在，我们可以使用以下命令绘制 Intel 的 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))
以下图显示了前述命令行的输出：
预测
现在，我们可以合理地确信我们的风险模型工作正常，我们还可以生成 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%分位数（在qnorm(0.99)中输入）计算 99%VaR。因此，下一期的一个月 99%VaR 估算值为qnorm(0.99)*0.1168 = 0.2717。因此，99%的概率下，月度回报将在-27%以上。
概要
在本章中，我们使用 R 解决了时间序列分析中的选定问题。我们涵盖了表示时间序列数据的不同方式，使用 ARMA 模型预测房价，通过协整关系改进了基本的最小方差对冲比率，并采用 GARCH 模型进行风险管理。在下一章中，您将学习如何使用 R 构建一个最优投资组合。
第二章：投资组合优化
到目前为止，我们已经熟悉了R语言的基础知识。我们知道如何分析数据、调用其内置函数，并将它们应用于时间序列分析中的选定问题。在本章中，我们将使用并扩展这些知识，讨论一个重要的实际应用：投资组合优化，换句话说，就是证券选择。本节内容涵盖了投资组合优化背后的理念：数学模型和理论解法。为了提高编程技能，我们将使用实际数据逐行实现一个算法来解决一个真实世界的例子。我们还将使用已编写好的 R 包来处理相同的数据集。
想象一下，我们生活在一个热带岛屿上，并且只有 100 美元可以投资。岛上的投资机会非常有限；我们只能将全部资金投资于冰淇淋或雨伞。根据天气变化，回报如下：
| 天气 | 冰淇淋 | 雨伞 |
|---|---|---|
| 晴天 | 120 | 90 |
| 阴天 | 90 | 120 |
假设天气下雨或晴天的概率相同。如果我们无法预见或改变天气，显然这两个选项是等价的，并且通过投资其中任何一个，我们将获得 5%的预期回报[(0.5×120+0.5×90)/100-1=0.05]。
如果我们可以在冰淇淋和雨伞之间分配我们的资金呢？那么我们应该分别投资 50 美元在这两个选项中。这个投资组合是无风险的，因为无论发生什么，我们通过一个资产赚取 45 美元，另一个赚取 60 美元。预期回报仍为 5%，但现在是有保障的，因为(45+60)/100-1=0.05。
投资组合优化的主要概念（该领域因哈里·马科维茨在 1990 年获得诺贝尔奖）可以通过这个例子来理解。根据投资产品之间的相关性，我们可以降低投资组合的风险（在此情况下通过方差来衡量），并且仍然获得预期的回报。
从数学上更精确地讲，设X和Y为具有有限方差的随机变量 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_01.jpg 和 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_02.jpg。它们的凸组合或仿射组合的方差表示如下的二次函数：
对于不同的相关性值，这个二次函数看起来像下图所示：
方差（作为风险的度量）只有在X和Y之间的相关性为-1 或+1，并且X和Y的方差不同的情况下，才能完全消除。否则，具有最优权重的投资组合的方差将（以绝对非平凡的方式）依赖于所有三个参数 (https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_06.jpg, https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_05.jpg, 和 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_07.jpg），正如我们在*定理（拉格朗日）*章节中将看到的那样。
均值-方差模型
马科维茨的均值-方差模型（*Markowitz, H.M. (1952 年 3 月)）*实际上是在更高维度中的冰淇淋/伞业务。对于数学公式化，我们需要一些定义。
它们的解释如下：
-
通过资产 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_08.jpg，我们指的是一个具有有限方差的随机变量。
-
通过投资组合，我们指的是资产的组合：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_09.jpg，其中https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_10.jpg，和https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_11.jpg。组合可以是仿射的或凸的。在仿射情况下，权重没有额外限制。在凸的情况下，所有权重都为非负数。
-
通过优化，我们指的是选择最佳的https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_12.jpg系数（权重）以使我们的投资组合满足我们的需求（即，在给定预期收益下具有最小风险，或者在给定风险水平下具有最高预期收益，等等）。
设https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_13.jpg为具有有限方差的随机收益变量，https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_14.jpg为它们的协方差矩阵，https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_15.jpg 和https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_16.jpg。
我们将重点关注以下优化问题：
很明显，https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_22.jpg 是投资组合的方差，https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_23.jpg 是预期收益。对于权重的和，我们有 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_24.jpg，这意味着我们希望投资 1 单位现金。（我们也可以考虑添加https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_25.jpg条件，这意味着不允许做空。）这些问题将在以下几点中详细解释：
-
第一个问题是找到一个最小风险的投资组合。如果没有无风险资产，这可能是一个非平凡的问题。
-
第二个问题是最大化给定方差水平下的预期收益。
-
一种略有不同的方法是找到在期望收益的指定水平下具有最小方差的投资组合。
-
第四个问题是最大化一个简单的效用函数 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_26.jpg，其中λ是风险容忍系数；它是一个任意数字，表示我们对风险的态度。它实际上与第一个问题相同。
-
在第五个问题中，Y是一个n+1 资产（例如，一个指数），我们不能购买或不想购买，但希望复制它。其他类似的问题可以用同样的方式来构建。
很明显，第二个问题是带有二次约束的线性优化问题；其他所有问题都是带有线性约束的二次函数。正如我们稍后所看到的，这是一个重要的区别，因为线性约束可以很容易地处理，而二次约束则更难处理。在接下来的两个章节中，我们将重点讨论这些问题的复杂性和可能的解决方案。
解的概念
在过去 50 年里，许多优秀的数值优化算法已经被开发出来，并且这些算法在二次函数的情况下特别有效。正如我们在前一节中看到的，我们只有二次函数和约束；因此，这些方法（在 R 中也有实现）可以在最坏情况下使用（如果没有更好的方法）。
然而，数值优化的详细讨论超出了本书的范围。幸运的是，在线性和二次函数及约束的特殊情况下，这些方法是不必要的；我们可以使用 18 世纪的拉格朗日定理。
定理（拉格朗日）
如果https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_27.jpg和https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_28.jpg，（其中https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_29.jpg）具有连续偏导数，并且https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_30.jpg是*f(x)*相对于https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_31.jpg约束的相对极值点，其中https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_32.jpg。
然后，存在系数https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_33.jpg，使得https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_34.jpg
换句话说，函数https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_35.jpg的所有偏导数都是 0（Bertsekas Dimitri P. (1999)）。
在我们的案例中，条件也是充分的。二次函数的偏导数是线性的，因此优化问题转化为解线性方程组的问题，这是高中水平的任务（与数值方法不同）。
让我们看看，如何利用这个来解决第三个问题：
可以证明，这个问题等价于以下线性方程组：
（协方差矩阵中增加了两行两列，因此我们也有条件来确定两个拉格朗日乘子。）我们可以预期这个系统有唯一解。
值得强调的是，通过拉格朗日定理得到的结果不再是一个优化问题。就像在一维情况下，最小化二次函数会转化为求导并解线性方程组，这在复杂性上是微不足道的。现在让我们看看如何处理回报最大化问题：
很容易看出，拉格朗日函数对λ的导数就是约束本身。
为了看清楚这一点，求拉格朗日函数L的导数：
因此，这导致了非线性方程，这更像是一门艺术而非科学。
处理实际数据
值得注意的是，投资组合优化已完全集成到各种 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 超集 (www.quandl.com/USER_1KR/1KT) 获取了一些 IT 股票价格数据，该服务是一个公共平台，提供便捷的途径访问大量的量化数据。虽然该 URL 指向一个可下载的逗号分隔值（CSV）文件 (www.quandl.com/api/v1/datasets/USER_1KR/1KT.csv)，可以保存到磁盘并通过 read.csv 导入 R，但通过前面 URL 中包含的密钥，有一种更直观的方式来完成此操作：
> 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"):
如果您没有使用认证令牌，前面的警告信息将会出现。请访问 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参数可选地指定我们感兴趣的时间段，并且设置为从现在起的最后 5 年。 -
请参阅
?Quandl获取更多选项；例如，type可用于导入已经存在于某些时间序列对象中的数据，而不是原始的data.frame。
在新创建的 IT 变量上运行 str 命令，显示了 R 对象的内部结构，该对象当前包含一个 Date 字段和五个资产的价格，格式为数字。
在将 IT 的价格（不包括第一个 Date 列）赋值给 assets 后，我们逐行运行前面的 minvariance 函数体。首先，我们通过将每个值（除第一个值外）与前一个值（tail 与 head）相除，并对每个商计算 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
接下来，我们开始构建拉格朗日定理指定的线性方程组的左侧：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_41.jpg，在这里我们结合了协方差矩阵（cov）、按列数（ncol）重复的 1（rep）和收益的均值（colMeans），并将其作为行（rbind）合并。
> 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
请注意，我们已将结果四舍五入到五位数字，以提高可读性。还要注意，微软（MSFT）和 AT&T 的平均回报是负数。现在，我们还将矩阵的最后两行（tail）作为新列（rbind）合并到左侧，以使其对于包含拉格朗日定理中额外零的线性系统变得完整（2x2 的matrix）：
> 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（在最小方差函数的参数中指定）；这是线性系统右侧向量的最后一个值。https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_42.jpg
> 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
请注意，除了微软和 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')
}
代码是一样的，唯一的不同是它会在最小和最大资产回报之间取一个数字（length = 100），计算不同回报值的最优投资组合的方差。然后我们可以绘制回报-方差对（s和rbase），以说明问题的解决方案。结果显示在以下截图中：
在方差-回报平面上，期望回报-最小方差曲线称为投资组合前沿。忽略其向下倾斜的部分（相同的方差可以通过更高的回报实现），我们得到有效前沿；没有理由选择有效前沿之外的投资组合。
众所周知，计算投资组合前沿只需对两个给定的回报水平进行计算，并将得到的投资组合结合起来，就能得到完整的前沿。
使用一些 R 包的内置函数，也可以实现类似的结果，而不需要太多编码。例如，fPortfolio包提供了一些有用的方法，已经可以直接应用于时间序列对象。为此，我们必须将原始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 中五只股票的回报率将如下图所示：
前面的前沿图可以通过绘制 portfolioFrontier 的结果进行交互式绘制：
> 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()
在前面的代码中，我们设置了一个特殊的 portfolioSpec S4 对象，并使用一个函数（solveRshortExact）来优化一个无限制的卖空组合。计算的结果（portfolioFrontier）通过 frontierPlot 绘制，并使用橙色圆点（pch = 19）；一些较小的（cex = 0.25）蒙特卡洛模拟点也添加到了图表中，并且在背景中有一个网格，如下图所示：
切点组合和资本市场线
当一个无风险资产https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_46.jpg被添加到模型中时会发生什么？如果https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_47.jpg和 X 是任何一个风险资产组合，那么https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_48.jpg，显然，https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_49.jpg。这意味着这些投资组合在均值-标准差平面上形成一条直线。任何位于这条直线上的投资组合都可以通过投资于 R 和 X 来获得。显然，对于 X 的最佳选择是这条直线与有效前沿相切的点。这个切点被称为市场组合或切点组合，并且在该点处，风险资产有效前沿的切线被称为资本市场线（CML），它包含了所有资产的有效投资组合。我们在均值-方差模型中要解决的最后一个问题是如何确定市场组合（或者等价地，CML）。
我们可以轻松地修改方差最小化代码来实现这一点。首先，如果我们添加一个无风险资产，则会在协方差矩阵中添加一个完全为零的行和列（其中 n 是包括无风险资产在内的资产数量）：
> n <- 6; mu <- 0.005
> Q <- cbind(cov(return), rep(0, n - 1))
> Q <- rbind(Q, rep(0, n))
然后，无风险回报（假设 rf 为 0.0001）被添加到回报向量中：
> 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
协方差矩阵中的噪声
当我们优化一个投资组合时，我们并没有实际的协方差矩阵和预期回报向量（这两个是均值-方差模型的输入）；我们使用观察值来估算它们，因此 Q、r 以及模型的输出也是随机变量。
不深入细节，我们可以说，这会导致模型中出现令人惊讶的巨大不确定性。尽管大数法则强力存在，最优投资组合权重有时会在https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_02_51.jpg之间变化。幸运的是，如果我们拥有几年的数据（每日回报），测量风险的相对误差仅为 20-25%。
当方差不足以满足需求时
方差作为风险度量虽然方便，但也存在一些缺点。例如，在使用方差时，回报的正向变化可能会被视为风险的增加。因此，已经开发出更复杂的风险度量方法。
例如，请参见以下关于各种方法的简短演示，这些方法应用于前面描述的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 表达式返回通过不同方法计算的投资组合权重，这些方法在本引言章节中没有讨论。有关详细信息，请参阅捆绑文档包中的文档，例如?portfolio，以及参考文献部分的相关文章和书籍章节。
概述
本章讨论了投资组合优化。在介绍了主要概念之后，我们介绍了马科维茨模型及其数学表达式。我们讨论了可能解决方案的方法，并实现了一个简单的算法，演示了这些方法如何在实际数据上运作。我们还使用了预编写的 R 包来解决同样的问题。我们广泛讨论了其他重要主题，如市场投资组合、协方差矩阵估计的不确定性，以及超越方差的风险度量。我们希望这是该主题的有用初步学习，并鼓励你进一步研究，或查阅下一章，该章涉及相关主题——资产定价模型。
第三章：资产定价模型
本章涉及了绝对定价问题（Cochrane 2005）以及如何根据资产的风险确定具有不确定支付的资产的价值。第二章，投资组合优化，在均值方差框架下通过分析资产回报建模了个人投资者的决策。本章则侧重于金融市场中是否可以存在均衡，所需的条件是什么，以及如何对其进行表征。将介绍两种主要方法——资本资产定价模型和套利定价理论，它们使用完全不同的假设和推理，但对回报演变给出了相似的描述。
根据相对定价的概念，基础产品的风险已经包含在其价格中，因此在衍生工具定价中不再发挥任何作用；这一点将在第六章，衍生品定价中呈现。无套利论证将在那里强制衍生资产与基础资产价格之间的一致性。
本章的目标是呈现资产回报与风险因子之间的关系。我们将解释如何从多个来源下载并清理数据。线性回归用于衡量依赖性，相关的假设检验则展示了结果的显著性。通过二步回归过程检验单因子指数模型，并展示结果的金融解释。
资本资产定价模型
解释资产价格的第一类模型使用了经济学考虑。利用前一章中提出的投资组合选择结果，资本资产定价模型（CAPM）回答了通过聚合理性投资者的决策，市场可以得出什么结论，以及均衡将如何演化的问题。Sharpe（1964）和Lintner（1965）证明了在以下假设条件下均衡的存在：
-
个人投资者是价格接受者
-
单期投资视野
-
投资仅限于可交易的金融资产
-
无税收且无交易成本
-
信息是免费的并且对所有投资者可得
-
投资者是理性均值-方差优化者
-
同质预期
在这些假设成立的世界里，所有投资者将持有相同的风险资产组合，即市场组合。市场组合包含所有证券，每个证券的比例是其市场价值占总市场价值的百分比。市场的风险溢价取决于所有市场参与者的平均风险厌恶程度。由此产生的均衡最著名的结果是市场风险溢价与单个证券风险之间的线性关系：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_A.jpg
E(r[i]) 是某证券的预期回报，r[f] 是无风险回报，E(r[m]) 是市场组合的预期回报。CAPM 中的风险由贝塔 β[i] 来衡量，贝塔是个体证券与市场的协方差与市场回报方差的函数：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_B.jpg
Cov[i,m] 是给定证券回报与市场回报之间的协方差，而 Var[m] 是市场回报的方差。
贝塔有多种解释。一方面，贝塔表示股票回报对市场组合回报的敏感度；另一方面，某个证券的贝塔表示该证券对市场组合所增加的风险。CAPM 模型指出，只有在系统性风险更高的情况下，市场才会提供更高的回报，因为非系统性风险可以通过多元化来分散，因此在此之后不能再支付风险溢价。
如果我们重新排列方程**(1)，就会得到所谓的证券市场线**（SML）的线性方程：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_C.jpg
CAPM 模型指出，在均衡状态下，所有证券都应位于 SML 线上；因此，尽管这些证券或投资组合可能不是有效的，这个方程对每个证券或投资组合依然成立。如果这个方程不成立，市场上就缺乏均衡。例如，如果某个证券的市场回报高于 CAPM 模型所预期的回报，那么每个投资者都需要调整自己投资组合的组成，以降低该证券的回报并使上述方程成立。
套利定价理论
套利定价理论 (APT)由罗斯（1977）提出，亦被用于金融领域来确定不同证券的回报。APT 模型指出，在均衡状态下，市场上不能存在套利机会，并且，资产的预期回报是多个随机因子的线性组合（Wilmott 2007）。这些因子可以是各种宏观经济因素或市场指数。在这个模型中，每个因子都有一个特定的贝塔系数：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_D.jpg
α[i] 是表示证券 i 的常数；β[ij] 是证券 i 对因子 j 的敏感度；F[j] 是系统性因子；而 e[i] 是证券的非系统性风险，均值为零。
APT 的一个核心概念是因子投资组合。因子投资组合是一个良好分散的投资组合，仅对一个因子做出反应，因此对所有其他因子的贝塔值为零，对该指定因子的贝塔值为 1。假设存在因子投资组合，可以通过套利论证来显示，任何良好分散的投资组合的风险溢价等于因子投资组合的风险溢价的加权和（Medvegyev-Száz 2010）。如果这个关系对每个良好分散的投资组合成立，则个别证券的预期回报将由因子的风险溢价（RP[j]）和其对该因子的敏感性（β[ij]）构成：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_E.jpg
如果 APT 模型中只有一个因子，即市场投资组合的回报，我们称该模型为指数模型。此外，如果α[i]为零，则我们将得到 CAPM 的精确定价公式。
CAPM 与 APT 的区别如下：
-
CAPM 是一个基于经济学考虑的均衡模型，而 APT 是一个基于套利论证的统计模型。
-
在 APT 模型中，如果拥有一个良好分散的投资组合，可以给出预期回报与贝塔值的关系，因此，通过在投资组合中包含大量资产，可以在实践中构建这种关系。而在 CAPM 模型中，所谓的市场投资组合是无法构建的。
-
CAPM 认为预期回报与贝塔值的关系适用于每个证券，而 APT 认为这种关系几乎适用于每个证券。
-
当市场存在错误定价时，在 APT 模型中，只要少数投资者改变投资组合结构，就足以得到证券的公平价格；而在 CAPM 模型中，每个投资者都必须如此操作。
贝塔估计
证券对某个因子的敏感性可以通过过去的价格波动来估计。我们将从单因子指数模型中估计贝塔值。首先，我们展示从不同来源收集和同步数据的过程，然后介绍简单的贝塔估计方法，最后建立一个线性回归模型。
数据选择
我们从 Quandl 下载某只股票（例如谷歌）和市场指数（标准普尔 500 指数）价格的时间序列，时间范围是 2009 年 6 月 1 日到 2013 年 6 月 1 日，如第二章所讨论：
> library(Quandl)
> Quandl.auth("yourauthenticationtoken")
> G <- Quandl('GOOG/NASDAQ_GOOG',
+ start_date = '2009-06-01', end_date = '2013-06-01')
得到的G是一个包含 6 个变量的变量，我们只需要收盘价的值：
> str(G)
'data.frame': 1018 obs. of 6 variables:
$ Date : Date, format: "2009-06-01" "2009-06-02" ...
$ Open : num 419 426 426 435 445 ...
$ High : num 430 430 432 441 447 ...
$ Low : num 419 423 424 434 439 ...
$ Close : num 427 428 432 440 444 ...
$ Volume: num 3323431 2626012 3535593 3639434 3681002 ...
> G <- G$Close
对标准普尔 500 指数的数据运行相同的代码，尽管现在我们处理的是调整收盘价：
> SP500 <- Quandl('YAHOO/INDEX_GSPC',
+ start_date = '2009-06-01', end_date = '2013-06-01')
> SP500 <- SP500$'Adjusted Close'
调整后的收盘价已经考虑了股息和拆股的影响。由于 Google 在此期间没有支付股息，也没有拆股，因此在本例中不需要进行此类调整。我们还需要无风险收益的时间序列，即 1 个月的美元 LIBOR 利率。虽然我们将使用每日回报，但 1 个月的利率可以视为短期利率，并且不容易受到隔夜利率中的随机噪音的影响。
> LIBOR <- Quandl('FED/RILSPDEPM01_N_B',
+ start_date = '2009-06-01', end_date = '2013-06-01')
> LIBOR <- LIBOR$Value
如你从之前的 Quandl 调用中看到的，每次数据都是从不同的数据提供商那里获取的。这也导致了数据结构的一些差异，因为我们有 Google 的Close值，S&P 500 的Adjusted Close值，以及LIBOR数据的Values值。向量的长度似乎也不相等：
> sapply(list(G, SP500, LIBOR), length)
[1] 1018 1008 1024
这意味着某些时间序列还包括其他时间序列中省略的日期。我们可以定义intersect函数来找到日期的交集，并在重新下载数据后仅筛选出这些数据单元格：
> G <- Quandl('GOOG/NASDAQ_GOOG',
+ start_date = '2009-06-01', end_date = '2013-06-01')
> SP500 <- Quandl('YAHOO/INDEX_GSPC',
+ start_date = '2009-06-01', end_date = '2013-06-01')
> LIBOR <- Quandl('FED/RILSPDEPM01_N_B',
+ start_date = '2009-06-01', end_date = '2013-06-01')
由于intersect函数只能应用于两个向量，我们调用Reduce函数来识别三个时间序列中共同的日期：
> cdates <- Reduce(intersect, list(G$Date, SP500$Date,LIBOR$Date))
现在，让我们简单地过滤所有三个数据框，只保留相关的单元格，以获得向量：
G <- G[G$Date %in% cdates, 'Close']
SP500 <- SP500[SP500$Date %in% cdates, 'Adjusted Close']
LIBOR <- LIBOR[LIBOR$Date %in% cdates, 'Value']
下载并清理数据后，你需要使用以下公式计算股票和市场指数的对数回报（r[t]）：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_F.jpg
S[t] 是第t天的市场价格。在 R 中，这可以表达为一个函数（请参见第二章，组合优化，了解更多细节）：
> logreturn <- function(x) log(tail(x, -1) / head(x, -1))
下一步，应该通过减去无风险的每日对数回报（r[ft]）来确定风险溢价。由于LIBOR利率是以货币市场为基础报价的——实际/360 天计数规则——且时间序列中包含的是百分比格式的利率，因此应使用以下公式：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_G.jpg
t 和 t-1 指的是日期，因此它们之间的差异是两个收盘值之间的天数，通常是 1 天，在我们的案例中，或者如果中间有非工作日，则会更多。使用以下命令，可以在 R 中轻松计算结果：
> rft <- log(1 + head(LIBOR, -1)/36000 * diff(cdates))
> str(rft)
num [1:1001] 1.81e-05 1.81e-05 1.81e-05 1.81e-05 5.42e-05 ...
我们通过计算共同日期之间的diff，已经计算出了(t+1)和t之间的差异。风险溢价（R[it]）可以表示为：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_H.jpg
简单的贝塔估算
一旦我们得到了两个时间序列；分别是个别资产（以 Google 为例）和市场（S&P 500）的风险溢价，就可以根据公式**(2)**计算贝塔值：
> cov(logreturn(G) - rft, logreturn(SP500) - rft) /
+ var(logreturn(SP500) - rft)
[1] 0.8997941
这也可以通过添加一个新函数来简化，用来描述风险溢价：
> riskpremium <- function(x) logreturn(x) - rft
> cov(riskpremium(G), riskpremium(SP500)) / var(riskpremium(SP500))
[1] 0.8997941
这种计算 beta 的方式与方程**(2)**有所不同，因为我们使用了风险溢价而不是回报。由于 CAPM 和 APT 都是单期模型，因此在两边都使用无风险回报进行修正不会影响结果。另一方面，在从时间序列估计 beta 时，我们必须决定是否使用回报或风险溢价，因为模型中的参数会有所不同，除非在常数无风险回报的情况下（Medvegyev-Száz 2010）。我们遵循前面描述的方法，就像遵循金融文献一样，但我们必须补充一点，Merryl Lynch 是通过回报计算 beta 的。
线性回归中的 Beta 估计
我们可以使用线性回归来估计 beta，其中解释变量是市场风险溢价（MRP），而因变量将是证券的风险溢价。因此，回归方程具有以下形式，这是证券特征线（SCL）的公式：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_I.jpg
我们将使用普通最小二乘法（OLS）估计来确定方程**(8)的线性回归模型。特征线的截距是α，即市场因素无法解释的股票回报部分。函数的斜率（方程(8)**）显示了对市场因素的敏感度，用 beta 来衡量。
我们可以使用 R 中内置的lm命令轻松计算回归模型：
> (fit <- lm(riskpremium(G) ~ riskpremium(SP500)))
Call:
lm(formula = riskpremium(G) ~ riskpremium(SP500))
Coefficients:
(Intercept) riskpremium(SP500)
0.0002078 0.8997941
我们不仅保存了结果，还打印了它们，因为我们添加了额外的括号。在模型的帮助下，也可以轻松绘制出 Google 的特征线，图表上显示 Google 的风险溢价与市场风险溢价的关系。
> plot(riskpremium(SP500), riskpremium(G))
> abline(fit, col = 'red')
下图展示了结果。横轴是 MRP，纵轴是 Google 股票的风险溢价：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933_03_1.jpg
根据 CAPM，α等于零，因此我们假设α[i]为 0，然后解除此限制。我们可以通过在模型中传入-1来强制α为零：
> fit <- lm(riskpremium(G) ~ -1 + riskpremium(SP500))
回归模型在 R 中的结果总结如下：
> summary(fit)
Call:
lm(formula = riskpremium(G) ~ -1 + riskpremium(SP500))
Residuals:
Min 1Q Median 3Q Max
-0.089794 -0.005553 0.000166 0.005520 0.117087
Coefficients:
Estimate Std. Error t value Pr(>|t|)
riskpremium(SP500) 0.90048 0.03501 25.72 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0124 on 1000 degrees of freedom
Multiple R-squared: 0.3982, Adjusted R-squared: 0.3976
F-statistic: 661.6 on 1 and 1000 DF, p-value: < 2.2e-16
较高的F-statistic值表明模型具有解释能力，beta 是显著的，原假设——beta 为零——在任何显著性水平下都应被拒绝。这些结果与 CAPM 一致。
如果我们通过放宽零假设α来进行测试，我们可以看到截距与零的差异不显著。较高的p-value值表明我们无法在任何常见的（超过 90%）显著性水平下拒绝原假设：
> summary(lm(riskpremium(G) ~ riskpremium(SP500)))
Call:
lm(formula = riskpremium(G) ~ riskpremium(SP500))
Residuals:
Min 1Q Median 3Q Max
-0.089999 -0.005757 -0.000045 0.005307 0.116883
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0002078 0.0003924 0.529 0.597
riskpremium(SP500) 0.8997941 0.0350463 25.674 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.01241 on 999 degrees of freedom
Multiple R-squared: 0.3975, Adjusted R-squared: 0.3969
F-statistic: 659.2 on 1 and 999 DF, p-value: < 2.2e-16
我们可以通过联合图检查残差，如下图所示。
> par(mfrow = c(2, 2))
> plot(fit)
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933_03_2.jpg
您还可以在PerformanceAnalytics包中找到内置函数CAPM.alpha和CAPM.beta，用于计算给定资产的阿尔法和贝塔参数。所需的参数是资产与基准资产的回报序列以及无风险利率。
模型测试
关于贝塔-回报关系的首次测试使用了二阶段线性回归法（Lintner 1965）。第一次回归估计了如上所述的个别证券的安全特征线和贝塔系数。在第二次回归中，证券的风险溢价是因变量，而贝塔系数是解释变量。原假设假定截距为零，曲线的斜率为市场风险溢价，市场风险溢价估算为样本的平均值。该测试可以通过增加一个解释变量来扩展：个别方差。
数据收集
我们将使用 2003 年至 2007 年间美国市场危机前时期的样本来展示测试。由于日数据包含更多的短期效应，我们将在从日时间序列计算出的月度回报上应用该测试。因此，我们需要更多股票的日价格时间序列；让我们下载 2003 年至 2007 年间按字母顺序排列的前 100 只标准普尔 500 指数成分股的价格：
> symbols <- c("A", "AA", "AAPL", "ABC", "ABT", "ACE", "ACN", "ACT", "ADBE", "ADI", "ADM", "ADP", "ADSK", "AEE", "AEP", "AES","AET", "AFL", "AGN", "AIG", "AIV", "AIZ", "AKAM", "ALL", "ALTR", "ALXN", "AMAT", "AMD", "AMGN", "AMP", "AMT", "AMZN", "AN", "ANF", "AON", "APA", "APC", "APD", "APH", "APOL", "ARG", "ATI", "AVB", "AVP", "AVY", "AXP", "AZO", "BA", "BAC", "BAX", "BBBY", "BBT", "BBY", "BCR", "BDX", "BEAM", "BEN", "BF.B", "BHI", "BIIB", "BK", "BLK", "BLL", "BMC", "BMS", "BMY", "BRCM", "BRK.B", "BSX", "BTU", "BXP", "C", "CA", "CAG", "CAH", "CAM", "CAT", "CB", "CBG", "CBS", "CCE", "CCI", "CCL", "CELG", "CERN", "CF", "CHK", "CHRW", "CI", "CINF", "CL", "CLF", "CLX", "CMA", "CMCSA", "CME")
请注意，前面的列表仅包含 96 只股票的名称，因为有四只股票在参考时间区间内缺失了太多数据。
让我们通过tseries包从统一的数据库中下载这些数据集：
> library(tseries)
> res <- lapply(symbols, function(symbol)
+ get.hist.quote(symbol, quote = "AdjClose", quiet = TRUE,
+ start = as.Date('2003-01-01'), end = as.Date('2007-01-01')))
因此，我们调用get.hist.quote函数，对每个symbol从默认的（Yahoo!）提供商下载调整后的收盘数据，不包含关于进度的任何详细信息（quiet）。请注意，抓取过程可能需要一些时间，并将产生 96 个时间序列的列表。现在，让我们也更新SP500和LIBOR，以适应新的时间区间，并定义新的公共日期交集：
> LIBOR <- Quandl('FED/RILSPDEPM01_N_B',
+ start_date = '2003-01-01', end_date = '2007-01-01')
> SP500 <- Quandl('YAHOO/INDEX_GSPC',
+ start_date = '2003-01-01', end_date = '2007-01-01')
> cdates <- intersect(LIBOR$Date, SP500$Date)
如上所述，我们需要一个月度数据集，而不是下载的日数据；让我们选取每个月的第一个值。为此，我们需要将公共日期的列表保存为Date格式：
> d <- data.frame(date = as.Date(cdates, origin = '1970-01-01'))
> str(d)
'data.frame': 998 obs. of 1 variable:
$ date: Date, format: "2003-01-02" "2003-01-03" ...
接下来，我们需要在相同的数据框中添加月份的日期以及拼接后的年份和月份：
> d$day <- format(d$date, format = '%d')
> d$my <- format(d$date, format = '%Y-%m')
现在我们只需对每组my（代表同一年同一个月）中的day变量应用min函数，day代表月份中的日期：
> (fds <- with(d, tapply(day, my, min)))
2003-01 2003-02 2003-03 2003-04 2003-05 2003-06 2003-07 2003-08
"02" "03" "03" "01" "01" "02" "01" "01"
2003-09 2003-10 2003-11 2003-12 2004-01 2004-02 2004-03 2004-04
"02" "01" "03" "01" "02" "02" "01" "01"
2004-05 2004-06 2004-07 2004-08 2004-09 2004-10 2004-11 2004-12
"03" "01" "01" "02" "01" "01" "01" "01"
2005-01 2005-02 2005-03 2005-04 2005-05 2005-06 2005-07 2005-08
"03" "01" "01" "01" "02" "01" "01" "01"
2005-09 2005-10 2005-11 2005-12 2006-01 2006-02 2006-03 2006-04
"01" "03" "01" "01" "03" "01" "01" "03"
2006-05 2006-06 2006-07 2006-08 2006-09 2006-10 2006-11 2006-12
"01" "01" "03" "01" "01" "02" "01" "01"
我们需要再次将结果与日期合并：
> (fds <- as.Date(paste(row.names(fds), fds, sep = '-')))
[1] "2003-01-02" "2003-02-03" "2003-03-03" "2003-04-01" "2003-05-01"
[6] "2003-06-02" "2003-07-01" "2003-08-01" "2003-09-02" "2003-10-01"
[11] "2003-11-03" "2003-12-01" "2004-01-02" "2004-02-02" "2004-03-01"
[16] "2004-04-01" "2004-05-03" "2004-06-01" "2004-07-01" "2004-08-02"
[21] "2004-09-01" "2004-10-01" "2004-11-01" "2004-12-01" "2005-01-03"
[26] "2005-02-01" "2005-03-01" "2005-04-01" "2005-05-02" "2005-06-01"
[31] "2005-07-01" "2005-08-01" "2005-09-01" "2005-10-03" "2005-11-01"
[36] "2005-12-01" "2006-01-03" "2006-02-01" "2006-03-01" "2006-04-03"
[41] "2006-05-01" "2006-06-01" "2006-07-03" "2006-08-01" "2006-09-01"
[46] "2006-10-02" "2006-11-01" "2006-12-01"
然后再次过滤res数据框，只保留上面识别出的日期：
> res <- lapply(res, function(x) x[which(zoo::index(x) %in% fds)])
然后，在将列表与时间序列合并之后，转换为常规的data.frame格式并为列命名就非常简单：
> res <- do.call(merge, res)
> str(res)
'zoo' series from 2003-01-02 to 2006-12-01
Data: num [1:48, 1:96] 17.8 15.3 12.1 12.5 15 ...
Index: Date[1:48], format: "2003-01-02" "2003-02-03" ...
> res <- as.data.frame(res)
> names(res) <- symbols
这将导致一个 48 行 96 列的数据框。我们仍然需要按列计算每只下载的股票的回报，但为此，rft也应该根据每个月的第一个值进行更新：
> LIBOR <- LIBOR[LIBOR$Date %in% fds, 'Value']
> rft <- log(1 + head(LIBOR, -1)/36000 * as.numeric(diff(fds)))
> res <- apply(res, 2, riskpremium)
让我们也将标准普尔 500 指数的值过滤为月度数据集：
> SP500 <- SP500[SP500$Date %in% fds, 'Adjusted Close']
建模 SCL
利用股票收益的时间序列，我们可以计算每只证券的贝塔。因此，我们将得到一个风险溢价向量，它是样本数据的平均值，以及一个包含贝塔的向量。
要估计的第二个回归模型如下：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_J.jpg
通过计算每只股票的riskpremium参数，并将其强制转换为data.frame，可以用基本循环一次性计算每只证券的贝塔和收益的均值：
> res <- apply(res, 2, riskpremium)
> res <- as.data.frame(res)
> r <- t(sapply(symbols, function(symbol)
+ c(beta = lm(res[, symbol] ~
+ riskpremium(SP500))$coefficients[[2]],
+ mean = mean(res[, symbol]))
+ ))
> r <- as.data.frame(r)
因此，遍历所有符号，让我们绘制计算得到的贝塔列表和风险溢价的平均值，如下图所示：
> plot(r$beta, r$mean)
> abline(lm(r$mean ~ r$beta), col = 'red')
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933_03_3.jpg
该模型可以描述如下：
> summary(lm(r$mean ~ r$beta))
Call:
lm(formula = r$mean ~ r$beta)
Residuals:
Min 1Q Median 3Q Max
-0.024046 -0.008783 -0.003475 0.006485 0.039731
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.009084 0.002429 3.740 0.000325 ***
r$beta 0.005528 0.001678 3.295 0.001413 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.01383 on 89 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.1087, Adjusted R-squared: 0.09873
F-statistic: 10.86 on 1 and 89 DF, p-value: 0.001413
根据上述结果，截距为正，但与零的差异不显著。SML 的斜率为 0.5528%——按月计算——略低于预期，因为根据零假设，它应该是该期间市场风险溢价的平均值：0.69%。然而，这一差异在统计上也是不显著的。基于测试，不能拒绝贝塔回报关系。
测试个体方差的解释力
该测试可以进一步发展，涉及作为第二个解释变量的非系统性风险。某个证券的个体风险可以通过以下方式计算：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_K.jpg
所以，首先我们必须计算方差向量，然后得到个体方差向量。要估计的回归方程如下：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_03_L.jpg
到目前为止，我们更新了用于计算贝塔和均值的上述循环，在r中：
> r <- t(sapply(symbols, function(symbol) {
+ stock <- res[, symbol]
+ beta <- lm(stock ~ riskpremium(SP500))$coefficients[[2]]
+ c(
+ beta = beta,
+ mean = mean(stock, na.rm = TRUE),
+ risk = var(stock, na.rm = TRUE) - beta² * var(SP500))
+ }))
> r <- as.data.frame(r)
尽管这个循环与之前的几乎相同，但大部分内容都根据DRY（不要重复自己）原则进行了重写和重新格式化。因此，首先我们将symbol的值存储在stock中，并在返回结果之前计算beta，并将其与c合并。现在，我们还为mean和var函数添加了na.rm = TRUE参数，以在计算之前移除可能的缺失值。我们的模型现在如下所示：
> summary(lm(r$mean ~ r$beta + r$risk))
Call:
lm(formula = r$mean ~ r$beta + r$risk)
Residuals:
Min 1Q Median 3Q Max
-0.023228 -0.009175 -0.003657 0.006817 0.036262
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.400e-02 3.711e-03 3.772 0.000285 ***
r$beta -1.743e-03 4.677e-03 -0.373 0.710293
r$risk -9.956e-08 5.798e-08 -1.717 0.089266 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.01381 on 93 degrees of freedom
Multiple R-squared: 0.1451, Adjusted R-squared: 0.1267
F-statistic: 7.891 on 2 and 93 DF, p-value: 0.0006833
有趣的是，新的参数将贝塔的回归系数改为了负值。然而，另一方面，风险参数在 95%的显著性水平上被证明是无关紧要的。由于资本资产定价模型（CAPM）认为对于可分散风险不需要支付风险溢价，零假设假定β[2]为零。在这里，我们无法拒绝这一假设。
米勒和斯科尔斯（1972）解释了第一次资本资产定价模型（CAPM）测试的结果——α与零显著不同，且斜率远低于市场风险溢价的平均值——并给出了统计学原因。作为第二次回归（贝塔）的解释变量，它是通过第一次回归估算得出的，因此包含统计误差。这种估计偏差导致了观察到的显著截距和比预期更加平坦的证券市场线（SML）。这一结论可以通过模拟收益进行检验。关于模拟的更多细节可以在接下来的两章中找到。
总结
本章中，资产收益的系统性风险通过它们对市场方差的贡献——贝塔系数来衡量。我们使用线性回归来量化这种关系，并进行了假设检验，以确认资本资产定价模型的相关结论。
第四章：固定收益证券
在第三章《资产定价模型》中，我们重点讨论了建立风险与其贝塔值、金融工具价格和投资组合之间关系的模型。第一个模型是 CAPM，采用了均衡方法，而第二个模型 APT 则建立在无套利假设的基础上。
固定收益组合管理的一般目标是建立一个具有特定风险/回报特征的固定收益证券组合。换句话说，投资组合经理的目标是将资金分配到不同的固定收益证券中，以最大化投资组合的预期回报，同时遵循既定的投资目标。
该过程包括收益率曲线、提前偿还行为和证券违约的动态建模。使用的工具包括时间序列分析、随机过程和优化。
固定收益证券的风险包括信用风险、流动性风险和市场风险等。前两者可以通过选择具有预定违约风险的证券来处理，例如，选择具有最低信用评级并具备适当流动性特征的证券。固定收益证券的市场风险通常通过久期、修正久期、基准久期或因子久期来衡量。这些都是衡量固定收益证券面临的利率风险的指标。本章重点讨论固定收益证券的市场风险。
衡量固定收益证券的市场风险
给定收益率曲线，获得固定收益证券现值的一般公式为：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_04_01.jpg，其中 T 是证券到期的时间，CF[t] 是证券在时间 t 的现金流，y[t] 是在时间 t 收到的现金流的折现率。即使到期收益率保持不变，债券的市场价格会随着时间的推移趋向其面值。这种价格变化是预期之中的，因此不被视为风险。市场风险来自利率的变化，这会导致再投资风险和清算风险。再投资风险影响的是息票支付的再投资利率，而清算风险则影响债券的市场价格。
利率变化对市场价格的影响是通过检查债券价格作为到期收益率(y)的函数来衡量的：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_04_02.jpg。由于https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_04_03.jpg，由https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_04_04.jpg引起的价格变化百分比可表示为：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_04_05.jpg，即！衡量固定收益证券的市场风险的二阶近似为！衡量固定收益证券的市场风险。当收益率按周期表示时，复利久期(D)、修正久期(D*)，和凸度的定义如下：
债券的定价公式显示了债券到期收益率(y)与其价格*(P*)之间的明显反向关系。由于久期与债券到期收益率的变化与其价格变化的关系密切，它是衡量债券利率风险的最重要指标。久期是债券加权平均到期时间。
示例 – 在 R 中的实现
假设有一只 10 年期债券，面值为 1,000 美元，到期时支付本金，年利率 8%，按季度支付票息，并假设收益率曲线在 10%的水平上使用连续复利计算。
为了计算上述描述的指标，我们将使用GUIDE包，它提供了一个图形用户界面，用于各种金融计算器和定价金融衍生品的交互式图表。因此，在以下示例中，大多数参数将比其他章节以更直观的方式进行设置。
安装并加载包后，可以通过GUIDE函数启动主程序：
> install.packages('GUIDE')
> library(GUIDE)
> GUIDE()
这将加载主窗口，显示一个菜单，用于访问该包的 55 个函数，如以下截图所示：
这些函数也可以通过直接的 R 命令在顶部菜单栏旁边调用。
公允价值由bondprice快速给出，金额为 867.28 美元。priceyield函数展示了折现率与债券价格之间的反向关系。债券的久期由bonddur确定，如以下截图所示：
该函数可以设置为反映年度或半年度的票息支付，且折现率的频率可以变化。该函数还允许计算修正久期。相同债券的凸度通过bondconv计算，如以下截图所示：
请注意，当贴现率以连续复利收益率表示时，凸度由以下公式给出：
duryield和durcoupon函数可用于评估收益率增加如何影响久期，以及较大的票息如何影响债券的久期。
久期与到期之间的关系通过durmaturity函数显示，如下截图所示：
计算了债券的久期或凸度后，投资组合的久期或凸度可以通过将投资组合各个组成部分的久期或凸度加权平均来轻松计算。
其他软件包，如maRketSim和termstrc，也包含了能够计算债券或整个固定收益投资组合的久期、修正久期和凸度的功能。
固定收益投资组合的免疫化
当投资组合不受利率变化的影响时，该投资组合就被视为已免疫。久期提供了一个很好的利率敏感性度量，因此通常用于免疫化投资组合。由于使用久期假设收益曲线是平坦的，并且收益曲线发生的是小的平行移动，因此免疫化投资组合受这些假设的限制，不受影响意味着随着收益变化，投资组合的价值仅会轻微变化。
有两种不同类型的免疫化策略：净资产免疫化和目标日期免疫化。
净资产免疫化
固定收益投资组合经理通常会对收益曲线未来的变化有一定的看法。假设一位投资组合经理预期短期内利率会上升。由于这对投资组合将产生不利影响，投资组合经理可以决定通过进入远期协议或利率互换，将投资组合的久期调整为零。这些工具可以改变投资组合的久期，并帮助将投资组合的久期调整为零，而无需清算整个投资组合。
投资组合经理的另一个目标是根据投资组合基准的久期来设定投资组合的久期。如果他们对市场走势的预期得到验证，这有助于超越投资组合的基准表现。
银行通常更关心保护其股本价值不受市场价格变化的影响。这通常通过将股本的久期设置为所需的水平来实现。
目标日期免疫化
假设我们考虑一个投资者，他有一个特定的负债现金流。通过构建一个固定收益证券的资产组合来实现投资者投资组合的免疫化，该资产组合的久期等于负债的久期。这种目标日期免疫化确保将来支付义务将从投资组合的资产中得到满足。该过程可以通过例如 genPortfolio.bond 函数来实现。
专项免疫化
专项免疫化是一种特殊的目标日期免疫化方式，其中资产的现金流与每个负债的组成部分相匹配。可以通过使用零息债券为负债组件提供资金来实现这一目标。
可转换债券定价
可转换债券通常由信用评级较低且增长潜力较大的公司发行。这些公司可以通过赋予债券持有人将债券转换为一定数量的公司普通股的权利（但没有义务），来降低其利息成本。投资者在享有转换为股权的潜在上涨空间的同时，还能通过债券提供现金流获得下行保护。公司从可转换债券转换时，公司的杠杆效应会减少，而转股时的股本稀释则是这一过程的权衡。
这些特征表明可转换债券的行为有三个不同的阶段：价内可转换债券（转换价格 < 股票价格）表现得像股票，平价可转换债券（转换价格 = 股票价格）被视为股票和债务的组合，而价外可转换债券（转换价格 > 股票价格）则表现出债务特征。即使在 Black-Scholes-Merton 模型框架下，为可转换债券定价也可能很复杂，但基本原则是分别对债券和期权进行定价。
假设我们考虑一只 5 年期的可转换债券，面值为 100 美元，年息 5%，年付利息，并且在到期时可将面值转换为 4 股普通股票。假设所有期限的无风险利率为 5%，该债券的信用利差为 2%，基础股票的价格为 20 美元，股票的波动率为 20%，股息收益率为零。R 可用于为该可转换债券定价。首先，我们定义将用于以下示例的今天日期：
> today <- Sys.Date()
接下来，我们设置交易和结算日期，并在给定平坦收益率曲线的情况下计算折现曲线的值（基于 times 参数，即当前在 0 到 10 之间的序列，步长为 0.1）：
> params <- list(tradeDate = today - 2,
+ settleDate = today,
+ dt = 0.25)
> times <- seq(0, 10, 0.1)
> dividendYield <- DiscountCurve(params, list(flat = 10e-6), times)
> riskFreeRate <- DiscountCurve(params, list(flat = 0.05), times)
前述的股息收益率、无风险利率，以及以下固定基础资产的价格和波动性将传递给后续的 Black-Scholes 过程，后者将为该债券建立二项式定价引擎：
> process <- list(
+ underlying = 20,
+ divYield = dividendYield,
+ rff = riskFreeRate,
+ volatility = 0.2)
我们还应该指定转换比例，这决定了债券持有人如果决定将债券转换为股票时，将获得多少普通股股份。债券的面值和信用利差也在此处指定：
> bondparams <- list(
+ exercise = "eu",
+ faceAmount = 100,
+ redemption = 100,
+ creditSpread = 0.02,
+ conversionRatio = 4,
+ issueDate = as.Date(today + 2),
+ maturityDate = as.Date(today + 1825))
具有年付息的债券：
> dateparams <- list(
+ settlementDays = 3,
+ dayCounter = "ActualActual",
+ period = "Annual",
+ businessDayConvention = "Unadjusted")
并将上述指定的参数传递给 ConvertibleFixedCouponBond 函数：
> ConvertibleFixedCouponBond(bondparams, coupon = 0.05, process, dateparams)
Concise summary of valuation for ConvertibleFixedCouponBond
Net present value : 107.1013
clean price : 107.06
dirty price : 107.1
accrued coupon : 0.041096
yield : 0.033848
cash flows :
Date Amount
2014-06-21 4.9589
2015-06-21 5.0000
2016-06-21 5.0073
2017-06-21 4.9927
2018-06-21 5.0000
2018-06-21 100.0000
如果不考虑可转换特性，债券的价值大约为 92 美元，而包含额外特性的债券价值为 107.1 美元。现在，让我们检查一下当我们开始将标的股票的价格从 1 提高到 30 时，净现值的变化：
> res <- sapply(seq(1, 30, 1), function(s) {
+ process$underlying = s
+ ConvertibleFixedCouponBond(bondparams, coupon = 0.05, process, dateparams)$NPV
+ })
> plot(1:30, res, type = 'l',
+ xlab = 'Price of the underlying stocks',
+ ylab = 'Net Present Value')
下图显示了标的股票价格与计算出的可转换债券价值之间的关系：
概要
在本章中，我们使用 R 来衡量固定收益投资组合的利率风险。我们介绍了 GUIDE 包中的部分函数，并应用了 RQuantLib 包中的可转换债券定价函数。在下一章中，你将学习如何使用 R 来估算即期收益率曲线。
第五章：估算利率期限结构
在上一章中，我们讨论了利率水平和期限结构的变化如何影响固定收益证券的价格。现在，我们将重点讨论利率期限结构的估算，这是金融学中的一个基本概念。它几乎在所有金融决策中都是一个重要的输入。本章将介绍通过三次样条回归估算期限结构的方法，并演示如何使用termstrc包和govbonds数据集来估算利率期限结构。
利率期限结构及相关函数
一只t年期的零息债券，面值为 1 美元，是一种在到期时支付 1 美元的证券，也就是说，t年后支付。设 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_03.jpg 为其市场价值，也称为t年期的折现因子。该函数 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_04.jpg 称为折现函数。基于无套利假设，通常假定 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_05.jpg，https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_06.jpg 是单调递减的，并且 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_07.jpg。还通常假定 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_06.jpg 是二次连续可微的。
设 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_08.jpg 表示t年期零息债券的连续复利年收益率；它定义为：
函数 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_10.jpg 称为（零息）收益率曲线。
设 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_11.jpg 表示瞬时远期利率曲线，或简称远期利率曲线，其中：
这里 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_13.jpg 是指在一个假设的远期贷款协议中，双方约定的利率，其中一方承诺在t年后，以非常短的期限，按签订合同时固定的利率将一笔金额借给另一方。
折现函数、收益率曲线和远期利率曲线相互决定，是利率期限结构的可能表示方式。利率期限结构可能与这些函数中的任何一个或全部相关。
估算问题
我们无法直接观察期限结构，但我们可以观察那些价格依赖于期限结构的工具的市场价格，从而估计期限结构。一个关于期限结构的良好信息来源是政府债券市场，在这里通常有大量的流动性证券交易，其价格完全依赖于期限结构。
假设有 n 种债券交易，其总价格（或脏价格）用 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_14.jpg 表示。有 m 个日期，在这些日期，至少有一种债券的持有者会收到支付。这些支付分别在 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_15.jpg 年后到期，其中 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_16.jpg 。矩阵 C 包含债券的现金流。我们将债券价格建模为债券现金流现值与一个正态分布误差项的和：
这里 d 是包含折现因子的向量 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_19.jpg，而 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_20.jpg 是包含误差项的向量。债券的观察市场价格可能与现金流现值有所不同，原因有两个：可能存在观测市场价格的测量误差和/或存在轻微的市场不完美，比如交易成本，这允许理论价格与市场价格之间存在小的差异，而这种差异并不构成套利机会。误差项的方差可能因债券不同而有所差异：
这里，https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_22.jpg 是一个 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_23.jpg 正半定对角矩阵。合乎逻辑的假设是，债券价格误差项的标准差与其买卖价差成正比，也就是债券的买价与卖价之间的差异。因此，https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_24.jpg 通常被选为债券 i 的买卖价差的平方。
方程式*(1)*看起来像一个典型的线性回归，然而，它通常不能直接估计，因为观察值（债券价格）的数量通常少于需要估计的系数数量。由于这个原因，我们需要对期限结构进行建模，以减少需要估计的参数数量，并确保最终的期限结构估计是合理的。
通过线性回归估计期限结构
假设折现函数可以表示为一个线性组合，其中的 https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_25.jpg 函数是二次连续可微的函数，如下所示：
哪里
我们可以通过广义最小二乘法估计权重https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_28.jpg。我们稍后将讨论函数https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_29.jpg的选择。估计的折现函数是估计权重https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_30.jpg的函数。
设D表示一个https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_32.jpg矩阵，其元素https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_33.jpg是https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_34.jpg，而https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_35.jpg是包含权重https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_28.jpg的向量。因此https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_36.jpg并且
这是一个线性回归模型，约束条件为https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_38.jpg，其可以表示为如下：
在方程（2）中，权重的 GLS 估计在方程（3）约束下为
其中
三次样条回归
如果我们希望估计能够得到合理的折现函数，我们需要仔细选择函数https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_29.jpg。典型的折现函数是非线性的。它是一个单调递减的函数，并且在无限远处渐近于零。因此，拟合一条直线并不是一个好主意。我们可以尝试拟合一个高阶多项式来逼近折现函数，但这也不是一个令人满意的解决方案。如果我们拟合低阶多项式，它们通常不够灵活，拟合效果也不好，尤其是在短期到期的情况下。如果我们拟合高阶多项式，它们可能拟合得较好，但往往会在长期到期时产生剧烈波动，特别是在到期债券相对较少的情况下。这些剧烈波动通常会导致不现实的期限结构估计。
样条函数是帮助解决这个问题的函数，因为它们的灵活性可以在需要的地方局部增加，而无需提高估计函数的多项式阶数。通过将三次样条拟合到折现函数来估计期限结构的做法是由麦卡洛克于 1971 年首次提出的。
三次样条是实函数，其定义域是实数线的一个区间。定义域https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_47.jpg被所谓的节点点https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_43.jpg分为子区间，其中https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_48.jpg。三次样条函数在每个子区间内是三次多项式，这些三次多项式在节点点处连接，使得样条函数在https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_44.jpg上连续且二次连续可导。每个三次样条函数在https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_44.jpg上和给定的节点点https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_43.jpg可以表示为https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_45.jpg基样条函数的线性组合，这些基样条函数是同一节点点上的三次样条。因此，如果我们想为折扣函数拟合一个三次样条，我们只需选择https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_29.jpg作为三次样条基函数，这一点我们将在govbonds数据集中的德国政府债券数据中演示。
> data(govbonds)
> str(govbonds[['GERMANY']])
List of 8
$ ISIN : chr [1:52] "DE0001141414" "DE0001137131" "DE0001141422" "DE0001137149" ...
$ MATURITYDATE: Date[1:52], format: "2008-02-15" "2008-03-14" ...
$ ISSUEDATE : Date[1:52], format: "2002-08-14" "2006-03-08" ...
$ COUPONRATE : num [1:52] 0.0425 0.03 0.03 0.0325 0.0413 ...
$ PRICE : num [1:52] 100 99.9 99.8 99.8 100.1 ...
$ ACCRUED : num [1:52] 4.09 2.66 2.43 2.07 2.39 ...
$ CASHFLOWS :List of 3
..$ ISIN: chr [1:384] "DE0001141414" "DE0001137131" "DE0001141422" "DE0001137149" ...
..$ CF : num [1:384] 104 103 103 103 104 ...
..$ DATE: Date[1:384], format: "2008-02-15" "2008-03-14" ...
$ TODAY : Date[1:1], format: "2008-01-30"
该数据集包含 52 只德国债券的信息，我们将集中关注发行和到期日期、价格以及提供的现金流。要创建一个类似的数据集以便进一步分析，请参见?govbonds的示例。
首先，我们使用prepro_bond函数预处理债券数据集，该函数返回现金流、到期时间、到期收益率、基于久期的权重矩阵、脏价格和应计利息向量等值：
> prepro <- prepro_bond('GERMANY', govbonds)
一个重要的决策是设置节点点的数量并确定它们的位置。第一个和最后一个节点点分别为零和 T，其他节点点通常选择使得每个子区间大致有相同数量的债券到期。设置节点点的数量并不那么简单。它将决定需要估计的参数数量，并将显著影响估计的期限结构。可以通过将K设置为1开始估计过程，然后将其增加一个，重复估计直到拟合优度有显著改善，并且估计的期限结构表现良好。或者，可以遵循 McCulloch 提出的经验法则，节点点的数量大约为https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_05_46.jpg。我们现在将在以下命令中演示这一点，借助从prepro对象分解出的到期矩阵：
> m <- prepro$m[[1]]
让我们定义德国债券的数量为n，并将其平方根的四舍五入值定义为s（节点点数量），结果为7：
> n <- ncol(m)
> s <- round(sqrt(n))
如果s的结果小于三，确定理想的节点点就很容易了。我们将使用以下命令：
> c(floor(min(y[, 1])), max(m[, ncol(m)]))
在这里，我们已经识别出到期矩阵中第一列的最小值（min）和最后一列的最大值（max），并将这些结果四舍五入到下方的最大整数（floor）。
如果s大于三，则第一个和最后一个结点点的定义与前面的命令行一样，其他结点点则通过一些辅助向量计算，这些辅助向量的长度为s-3，如以下命令所示：
> i <- 2:(s-2)
> h <- trunc(((i - 1) * n) / (s - 2))
> theta <- ((i - 1) * n) / (s - 2) - h
i向量仅仅持有从2到5的序列，从中我们计算了到期矩阵的列（h）的索引，这些索引将用于搜索其他结点点。这里theta作为权重使用。
> apply(as.matrix(m[, h]), 2, max) ++ theta * (apply(as.matrix(m[, h + 1]), 2, max) – + apply(as.matrix(m[, h]), 2, max))
在这里，我们找到了到期矩阵中每一列的最大值，并在以下输出中加入了theta加权的h+1列和h列最大值的差异：
DE0001135101 DE0001141463 DE0001135218 DE0001135317
1.006027 2.380274 5.033425 9.234521
现在我们使用c函数将先前计算出的最小（min）值和最大（max）值与前面的结果连接起来，这些结果是从前面的代码块中复制的，用于识别所有结点点：
> c(floor(min(y[, 1])), apply(as.matrix(m[, h]), 2, max) + theta * (apply(as.matrix(m[, h + 1]), 2, max) - apply(as.matrix(m[, h]), 2, max)), max(m[, ncol(m)]))
DE0001135101 DE0001141463 DE0001135218 DE0001135317
0.0000 1.006027 2.380274 5.033425 9.234521 31.44657
应用的 R 函数
尽管我们已经在前面的例子中使用了一些termstrc包中的函数，演示了如何确定理想的结点点数量并指定它们，但这一过程也可以通过一些进一步的 R 函数更轻松地完成，如以下命令行所示：
> x <- estim_cs(govbonds, 'GERMANY')
> x$knotpoints[[1]]
DE0001135101 DE0001141463 DE0001135218 DE0001135317
0.0000 1.006027 2.380274 5.033425 9.234521 31.44657
首先，我们使用了estim_cs函数，该函数基于三次样条（Ferstl-Haydn, 2010）估算票息债券的期限结构，并返回一个包含knotpoints名称的结点点列表。请注意，estim_cs是与列表一起工作的——就像包中的大多数函数一样——这就是为什么x$knotpoints返回了一个列表，我们只检查了其中第一个元素，该元素与我们在前一部分手动计算的值相同。
前述函数返回了一些其他有用的值，默认情况下会导致以下命令行：
---------------------------------------------------
Estimated parameters and robust standard errors:
---------------------------------------------------
[1] "GERMANY:"
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
alpha 1 1.9320e-02 1.5230e-02 1.2686 0.2111
alpha 2 -8.4936e-05 3.7926e-03 -0.0224 0.9822
alpha 3 -3.2009e-04 1.1359e-03 -0.2818 0.7794
alpha 4 -3.7101e-04 3.9074e-04 -0.9495 0.3474
alpha 5 7.2921e-04 9.9560e-05 7.3243 3.375e-09 ***
alpha 6 2.0159e-03 1.3019e-04 15.4843 < 2.2e-16 ***
alpha 7 -4.1632e-02 4.5903e-03 -9.0696 1.011e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
拟合优度可以通过summary函数展示，就像其他 R 模型一样：
---------------------------------------------------
Goodness of fit:
---------------------------------------------------
GERMANY
RMSE-Prices 0.198573
AABSE-Prices 0.131036
RMSE-Yields (in %) 0.130108
AABSE-Yields (in %) 0.057223
零息收益率曲线及其置信区间可以通过简单地对x对象调用plot来轻松展示，其中的结点点也被包含在内。
前面的图表展示了估算的零息收益率曲线和数据库中各个债券的到期收益率。到期时间最短的两个债券是异常值，因为它们在到期前可能流动性较差。我们看到，估算的收益率曲线与 10 年的到期收益率非常接近。对于较长的到期时间，估算的零息收益率通常高于相应的票息债券到期收益率。这可能是拟合不完美的结果，或者可以通过这些债券不是零息债券而是票息债券来解释。
通过使用par设置mfrow, 我们还可以在同一框架中绘制两个不同的图形（我们还设置了multiple=TRUE，这样图形将在不等待用户输入的情况下渲染）。例如，我们可以使用以下命令绘制x对象中的折现曲线和远期曲线：
> par(mfrow = c(2,1))
> plot(x$discount, multiple = TRUE)
> plot(x$forward, multiple = TRUE)
在前面的图中，我们可以看到，估计的折现函数和远期利率曲线表现良好，不存在套利机会（折现函数单调递减，远期利率曲线不会产生不现实的值和波动）。
三次样条估计的期限结构通常能给出良好的估计。然而，有时估计的期限结构并不理想（远期利率曲线波动较大）。在这种情况下，可以使用非线性样条回归或简约收益率曲线模型，但这些内容超出了本章的范围。
通过estim_nss函数或YieldCurve包，您还可以使用包括 Nelson/Siegel、Diebold/Li、Svensson 和调整后的 Svensson 方法等更多资源。
总结
本章讨论了通过三次样条回归进行期限结构估计的方法，并演示了如何在 R 中估计利率的期限结构。在简要介绍期限结构和利率的理论后，还讨论了最基本的方法，如线性回归模型及相关问题，本章详细概述了三次样条回归模型在 R 中的实现，并提到了一些已经发布的 R 函数和包，以应对具有更复杂预期的任务。
第六章：衍生品定价
衍生品是其价值源自（或依赖于）另一产品价值的金融工具，称为标的资产。衍生品的三种基本类型是远期和期货合约、掉期和期权。在本章中，我们将重点讨论后者，并展示如何在 R 中处理基本的期权定价模型及一些相关问题。我们将从概述如何在 R 中使用连续的黑-斯科尔斯模型和二项式 Cox-Ross-Rubinstein 模型开始，然后继续讨论这些模型之间的联系。此外，通过计算和绘制希腊字母，我们将展示如何分析期权所涉及的最重要的市场风险。最后，我们将讨论隐含波动率的含义，并通过实际市场数据绘制波动率微笑来说明这一现象。
与期货或掉期相比，期权的最重要特点是你不能确定交易（买入或卖出标的资产）是否会发生。这一特性使得期权定价更加复杂，并要求所有模型对标的资产未来价格走势做出假设。我们在此讨论的两个模型在这些假设上有所不同：黑-斯科尔斯模型采用连续过程，而 Cox-Ross-Rubinstein 模型则采用离散随机过程。然而，其他假设非常相似，我们会看到它们的结果非常接近（此外，原则上是相同的）。
黑-斯科尔斯模型
黑-斯科尔斯模型的假设（Black 和 Sholes, 1973，参见 Merton, 1973）如下：
-
标的资产（S）的价格遵循几何布朗运动：https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_06_01.jpg
这里 µ（漂移）和 σ（波动率）是常数参数，W 是标准维纳过程。
-
市场是无套利的。
-
标的资产是没有分红的股票。
-
可以以任何（甚至是分数）数量买入和（做空）卖出标的资产。
-
没有交易成本。
-
短期利率（r）已知且在一段时间内保持不变。
该模型的主要结果是，在这些假设下，欧洲看涨期权（c）的价格具有封闭形式：
这里 X 是行权价，T-t 是期权的到期时间，N 表示标准正态分布的累积分布函数。给出期权价格的方程通常称为黑-斯科尔斯公式。从期权-看涨-平价关系可以容易地看出，具有相同参数的欧洲看跌期权（p）的价格由下式给出：
现在考虑 2013 年 6 月的谷歌股票看涨和看跌期权，期权到期日为 2013 年 9 月（即，距到期有 3 个月）。假设当前标的股票价格为 900 美元，行权价为 950 美元，谷歌的波动率为 22%，无风险利率为 2%。我们将使用fOptions包中的GBSOption函数计算看涨期权的价值。除了之前讨论的参数外，我们还必须设置持有成本（b）；在原始的 Black-Scholes 模型中（假设标的没有支付股息），它等于无风险利率。
> library(fOptions)
> GBSOption(TypeFlag = "c", S = 900, X =950, Time = 1/4, r = 0.02,
+ sigma = 0.22, b = 0.02)
Title:
Black Scholes Option Valuation
Call:
GBSOption(TypeFlag = "c", S = 900, X = 950, Time = 1/4, r = 0.02,
b = 0.02, sigma = 0.22)
Parameters:
Value:
TypeFlag c
S 900
X 950
Time 0.25
r 0.02
b 0.02
sigma 0.22
Option Price:
21.79275
Description:
Tue Jun 25 12:54:41 2013
这个延长的输出返回传递的参数，结果显示在Option Price标签下方。将TypeFlag设置为p将计算看跌期权的价格，现在我们只关心结果（可以在price槽中找到——有关更多细节，请查看对象的str）而不需要文本输出：
> GBSOption(TypeFlag = "p", S = 900, X =950, Time = 1/4, r = 0.02, sigma = 0.22, b = 0.02)@price
[1] 67.05461
和上一章一样，我们还可以选择使用GUIDE包提供的更易用的计算器来计算上述值。运行blackscholes()函数将触发一个模态窗口，其中有一个表单，我们可以在其中输入相同的参数。请注意，函数使用股息收益率代替持有成本，在本例中持有成本为零。
Cox-Ross-Rubinstein 模型
Cox-Ross-Rubinstein（CRR）模型（Cox, Ross 和 Rubinstein，1979）假设标的资产的价格遵循离散的二项过程。价格在每个周期可能上涨或下跌，因此根据二项树变化，如下图所示，其中u和d是固定的乘数，衡量价格上涨和下跌时的变化。CRR 模型的一个重要特性是u=1/d，且树是重组的；也就是说，如果价格先上涨再下跌或反之，两个周期后的价格将相同，如下图所示：
要构建二项树，我们首先必须决定我们要模拟多少步（n）；即期权到期时间将被划分成多少步。或者，我们可以确定树上每个时间步长的长度https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_06_a.jpgt，（以年为单位）：
如果我们知道标的资产的波动率（σ），则可以根据以下公式确定u和d参数：
因此：
在二项模型中定价期权时，我们需要确定标的资产的价格树，直到期权到期。然后，在获取所有可能的到期价格后，我们可以根据以下公式计算相应的期权价值：
为了使用二项模型确定期权价格，在每个节点，我们必须计算下一个两个可能期权值的期望值，然后进行折现。问题在于，折现时使用什么期望回报并不显而易见。窍门是，我们是通过假设的概率计算期望值，这使得我们可以使用无风险利率进行折现。这个概率叫做风险中性概率（p[n]），可以通过以下公式确定：
风险中性概率的解释是相当合理的：如果基础资产价格上涨的单期概率是 p[n]，那么基础资产的期望回报将是无风险利率。因此，使用 p[n] 计算的期望值可以按照 r 进行折现，期权在树中任何节点的价格都可以通过以下公式确定：
在上述公式中，g 是某一节点中期权的价格（它可以是看涨期权或看跌期权），g[u] 和 g[d] 是该衍生品在一个周期后两个可能节点的值。
为了在 R 中演示 CRR 模型，我们将使用与 Black-Scholes 公式中相同的参数。因此，S=900，X=950，σ=22% ，r=2%，b=2%，T-t=0.25。我们还需要设置 n，即二项树上的时间步数。为了演示，我们将使用一个 3 期模型：
> CRRBinomialTreeOption(TypeFlag = "ce", S = 900, X = 950,
+ Time = 1/4, r = 0.02, b = 0.02, sigma = 0.22, n = 3)@price
[1] 20.33618
> CRRBinomialTreeOption(TypeFlag = "pe", S = 900, X = 950,
+ Time = 1/4, r = 0.02, b = 0.02, sigma = 0.22, n = 3)@price
[1] 65.59803
值得注意的是，从二项模型得到的期权价格与之前计算的 Black-Scholes 价格接近（但并不完全相同）。除了最终结果，即期权的当前价格外，我们可能还会对整个期权树感兴趣：
> CRRTree <- BinomialTreeOption(TypeFlag = "ce", S = 900, X = 950,
+ Time = 1/4, r = 0.02, b = 0.02, sigma = 0.22, n = 3)
> BinomialTreePlot(CRRTree, dy = 1, xlab = "Time steps",
+ ylab = "Number of up steps", xlim = c(0,4))
> title(main = "Call Option Tree")
这里我们首先通过 BinomialTreeOption 计算了一个矩阵，使用给定的参数，并将结果保存在 CRRTree 中，然后将其传递给绘图函数，指定了 x 轴和 y 轴的标签，x 轴的范围设置为从 0 到 4，如下图所示。y 轴（上行步数）显示了基础资产总共上涨了多少次。下行步数定义为负的上行步数。
欧式看跌期权可以通过将之前代码中的 TypeFlag 更改为 pe 来类似显示：
https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/Image1.jpg
两个模型之间的联系
在应用了两种基本的期权定价模型后，我们给出了它们的一些理论背景。我们并不打算给出详细的数学推导，而是旨在强调（并随后在 R 中加以说明）这两种方法的相似性。连续期权定价和二项式期权定价背后的金融思想是相同的：如果我们通过持有适当数量的标的资产来完美对冲期权，这意味着我们创建了一个无风险投资组合。由于市场应该是无套利的，因此无风险投资组合的收益必须等于无风险利率。一个重要的观察是，正确的对冲比率是每个期权持有https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_06_16.jpg标的资产。因此，这个比率是期权价值相对于标的资产价格的偏导数（或其在二项模型中的离散对应物）。这个偏导数被称为期权的 delta。两种模型之间的另一个有趣的联系是，delta 对冲策略和相关的无套利论证得出了相同的定价原理：衍生品的价值是其未来可能价值的风险中性期望值，按无风险利率贴现。这个原理在二项树模型中非常容易处理，我们逐节点计算贴现的期望值；然而，连续模型也有相同的逻辑，即使期望值的数学计算更为复杂。这就是为什么我们只给出了这个论证的最终结果，即布莱克-斯科尔斯公式。
现在我们知道这两种模型具有相同的定价原则和思想（delta 对冲和风险中性估值），但我们也观察到它们的数值结果并不相等。原因是描述标的资产价格波动的随机过程并不相同。尽管如此，它们非常相似；如果我们像在*《Cox-Ross-Rubinstein 模型》一节中那样，根据波动率参数来确定u和d*的值，那么二项过程就可以逼近几何布朗运动。因此，随着时间步数的增加（或等效地，减少步长），二项模型的期权价格会趋近于布莱克-斯科尔斯模型的期权价格。
为了说明这种关系，我们将计算二项模型中期权价格随时间步数增加的变化。在下图中，我们将结果与期权的布莱克-斯科尔斯价格进行比较：
该图是通过一个循环生成的，该循环将N从1运行到200，以计算带有固定参数的CRRBinomialTreeOption：
> prices <- sapply(1:200, function(n) {
+ CRRBinomialTreeOption(TypeFlag = "ce", S = 900, X = 950,
+ Time = 1/4, r = 0.02, b = 0.02, sigma = 0.22, n = n)@price
+ })
现在prices变量包含了200个计算值：
> str(prices)
num [1:200] 26.9 24.9 20.3 23.9 20.4...
让我们还用广义布莱克-斯科尔斯期权计算该期权：
> price <- GBSOption(TypeFlag = "c", S = 900, X = 950, Time = 1/4, r = 0.02, sigma = 0.22, b = 0.02)@price
并在一个联合图中显示这些价格，GBS 选项以红色渲染：
> plot(1:200, prices, type='l', xlab = 'Number of steps',
+ ylab = 'Prices')
> abline(h = price, col ='red')
> legend("bottomright", legend = c('CRR-price', 'BS-price'),
+ col = c('black', 'red'), pch = 19)
希腊字母
理解期权可能涉及的风险类型对所有市场参与者至关重要。希腊字母的理念是衡量不同类型的风险；它们表示期权对不同因素的敏感性。普通看涨期权的希腊字母有：delta (https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_06_a.jpg，对基础资产价格的敏感性)，gamma (https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_06_b.jpg，delta 对基础资产价格的敏感性，delta 的 delta)，theta (https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_06_c.jpg，对时间的敏感性)，rho (https://github.com/OpenDocCN/freelearn-quant-zh/raw/master/docs/intd-r-quant-fin/img/0933OS_06_d.jpg，对无风险利率的敏感性)，以及vega (V，对波动率的敏感性)。从数学角度来看，所有希腊字母都是衍生品价格的偏导数：
可以通过GBSGreeks函数轻松计算每个期权的希腊字母：
> sapply(c('delta', 'gamma', 'vega', 'theta', 'rho'), function(greek)
+ GBSGreeks(Selection = greek, TypeFlag = "c", S = 900, X = 950,
+ Time = 1/4, r = 0.02, b = 0.02, sigma = 0.22)
+ )
delta gamma vega theta rho
0.347874404 0.003733069 166.308230868 -79.001505841 72.82355323
分析在市场参数变化时某个希腊字母如何变化通常是有用的。这种分析可能帮助我们更好地理解风险。例如，看涨期权的 delta 作为基础资产价格的函数是一个 S 形的递增曲线，范围从 0 到 1。这些特性始终有效，但如果时间流逝，我们接近期权到期日，曲线会变得越来越陡峭（参见下图）。其含义如下：如果看涨期权被行使的可能性很大，那么它就非常类似于一个长期远期合约，因此，delta 接近 1。如果行使的机会非常低，持有看涨期权就像持有空头仓位，delta 接近 0。随着时间的推移，基础资产价格区间内行使期权的不确定性（即既不太可能，也不太不可能）逐渐缩小；因此，delta 的曲线变得越来越陡峭。为了说明这一行为，我们将绘制看涨期权的 delta，作为基础资产价格的函数，展示三种不同的到期时间。
为了计算 delta，我们运行两个循环：一个循环使用三个不同的时间值，另一个循环让S从500跑到1500：
> deltas <- sapply(c(1/4, 1/20, 1/50), function(t)
+ sapply(500:1500, function(S)
+ GBSGreeks(Selection = 'delta', TypeFlag = "c",
+ S = S, X = 950, Time = t, r = 0.02, b = 0.02, sigma = 0.22)))
结果的deltas包含 1001 行（对应S值）和三列（对应指定的时间），我们将在联合图中展示：
> plot(500:1500, deltas[, 1], ylab = 'Delta of call option',
+ xlab = "Price of the underlying (S)", type = 'l')
> lines(500:1500, deltas[, 2], col='blue')
> lines(500:1500, deltas[, 3], col='red')
> legend("bottomright", legend = c('t=1/4', 't=1/20', 't=1/50'),
+ col = c('black', 'blue', 'red'), pch = 19)
下图展示了看涨期权在三种不同到期时间下的 delta：
确定或绘制复杂期权策略的希腊字母非常相似。例如，计算跨式头寸（一个由相同参数的看涨期权和看跌期权组成的投资组合）的 delta，只需分别计算看涨期权和看跌期权的 delta，然后将它们相加。我们将绘制跨式头寸的delta，作为基础资产价格的函数。我们可能会观察到，图形与之前看涨期权的 delta 非常相似，但现在的 S 曲线范围从-1 到 1：
> straddles <- sapply(c('c', 'p'), function(type)
+ sapply(500:1500, function(S)
+ GBSGreeks(Selection = 'delta', TypeFlag = type, S = S, X = 950, Time = 1/4, r = 0.02, b = 0.02, sigma = 0.22)))
所以我们调用一个嵌套循环，让S从500到1500，对看涨和看跌期权都进行计算，保持其他参数不变，并将结果的德尔塔值保存在一个矩阵中。接下来的命令将这些行（看涨和看跌期权）的和呈现出来：
> plot(500:1500, rowSums(straddles), type='l',
+ xlab='Price of the underlying (S)', ylab = 'Delta of straddle')
结果图展示了一个跨式期权头寸的德尔塔与标的资产价格的关系，如下图所示：
隐含波动率
Black-Scholes 模型常常因为一些缺点而受到批评。一个重要的问题是该模型假设标的资产的波动率是常数，而现实中这一假设并不成立。此外，由于波动率无法直接观测到，volatility是模型中最难以校准的参数。由于这一困难，Black-Scholes 公式通常以间接的方式来估算volatility参数；我们观察期权的市场价格，然后根据其他已知参数，寻找使得 Black-Scholes 定价与市场价格相等的σ值。这个σ参数被称为期权的隐含波动率。正如 Riccardo Rebonato 所言，隐含波动率是“在错误的公式中放入错误的数字以获得正确的价格”（Rebonato, 1999, p.78）。
我们将通过一些 Google 期权来说明隐含波动率的计算。这些期权是到期日为 2013 年 9 月 21 日的看涨期权，行权价格从 700 美元到 1150 美元不等（共有 76 个不同的期权）。我们在 2013 年 6 月 25 日从finance.google.com收集了这些期权的卖出价格，并将它们放入 CSV 文件中。进行计算时，我们需要知道当天 Google 的股价为 866.2 美元。由于距离到期日还有 88 天，我们将使用 88/360 年作为Time参数。无风险利率和持仓成本假设保持为 2%。
首先，从 CSV 文件中加载 Google 的期权数据：
> goog <- read.csv('goog_calls.csv')
然后，针对数据集的每一行运行一个循环，使用给定的参数计算波动率：
> volatilites <- sapply(seq_along(goog$Strike), function(i)
+ GBSVolatility(price = goog$Ask.Price[i], TypeFlag = "c",
+ S = 866.2, X = goog$Strike[i], Time = 88/360, r = 0.02, b = 0.02))
volatilities变量是一个向量，包含了计算出来的值：
> str(volatilites)
num [1:76] 0.258 0.253 0.269 0.267 0.257...
这可以与行权价格进行比较：
> plot(x = goog$Strike, volatilites, type = 'p',
+ ylab = 'Implied volatiltiy', xlab = 'Strike price (X)')
因此，下图展示了不同执行价格的隐含波动率：
值得注意的是，针对 Google 期权计算的隐含波动率会根据行权价格的不同而变化。这与 Black-Scholes 模型的假设相悖，后者认为波动率是常数。观察到的隐含波动率模式（中等行权价的波动率较低）并不是唯一的，且在金融市场中相当常见。由于曲线的特定形态，这一现象被称为“波动率微笑”。
总结
在本章中，我们使用 R 语言通过 Black-Scholes 模型和 Cox-Ross-Rubinstein 模型来定价普通的香草期权。此外，我们还研究了这些期权的基本希腊字母以及隐含波动率。有关这些主题的金融背景的更多细节，请参见(Hull, 2011)。除了了解 fOptions 包中的一些工具外，我们还编写了几个用于模拟目的的循环和自定义函数。下一章将集中讨论如何通过各种模型管理信用风险，例如通过蒙特卡洛模拟和信用评分方法选择最佳信用组合。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
