《Mastering R for Quantitative Finance》Chapter1笔记

嗯，本来是要学习《Mastering R for Quantitative Finance》第一章的，但是作者建议先学习《Introduction to R for Quantitative Finance》（一下简称introduction）的第一章，那就从这里入手吧。

Introduction的第一章讲的是时间序列分析，主要包括：时间序列获取和处理、线性时间序列分析及其在房价预测

上的应用、使用协整考虑long-run trend来优化最小方差对冲比率（basic minimal variance hedge ratio），最后还介

绍如何使用波动率模型实现风险管理目标。以下开始内容部分笔记：

1、数据获取和处理

R中存放时间序列的数据类型可以是向量、矩阵、dataframe以及特定的ts，但是这些都有局限性，一个比较好的

数据类型是 zoo ,  xts ,或者 timeSeries 对象。这些数据类型或者对象的操作可能更加适合时间序列分析，比如做季

节模型，估计像差分这种操作，在以上三种对象中可能更容易首先一些吧。这三种对象，分别安装同名package即可。

library(quantmod)                 #加载quantmod包，下载股票价格用到的package
options(stringsAsFactors = FALSE) #选项设置，避免数据读取的时候把string转换成factor
getSymbols('AAPL',src = "yahoo",from="2002-01-01")
class(AAPL)                       #"xts" "zoo",xts是zoo的一种扩展
head(AAPL);tail(AAPL)
df=AAPL$AAPL.Adjusted             #选取其中的调整收盘价作为分析对象
plot(df,main='苹果股票价格趋势图',ylab='价格（美元）',xlab='日期')
df[which.max(df)];df[which.min(df)]##查看最高股价和最低股价

通常，研究股票我们最终关注的不是股票价格的绝对数值，而是收益率。

收益率往往研究两种：简单的收益率和连续复利收益率。简单收益率计算公式为：

(P(t)-P(t-1))/P(t-1)

连续复利收益率计算公式为：

log(P(t))-log(P(t-1))

R 中的计算方式是：

rtn_simp=diff(df)/lag(df,k=-1)   #简单收益率
rtn_cont=diff(log(df))   #连续复利收益率
summary(rtn_simp);summary(coredata(rtn_simp))  #习惯性summary一下
hist(rtn_simp*100, breaks=100, main = "简单收益率分布直方图",xlab="百分之");
hist(rtn_simp*100,freq=F, breaks=100, main = "简单收益率分布直方图",xlab="%"); # freq=F

我们爬取的数据时间窗口比较长，如果想要取其子集，可以使用window函数：

rtn_simp_2016 = window(rtn_simp, start = '2016-01-01', end = '2016-12-31')

2、线性时间序列分析和预测

线性时间序列分析是Box和Jenkins在1976年提出的方法，在数学和统计学的基础理论中，还是非常年轻的。主要

的方法就是ARIMA模型。模型建立主要遵从如下步骤：模型识别->模型估计->模型诊断。模型识别主要是定阶(order)，

即模型中滞后项的order，通常是通过作图或信息量准则进行定阶。模型估计就是通常的最小二乘或者最大似然估计。

模型诊断指的是残差项的白噪声检验。下面直接用数据分析来展示这个一般流程。

2.1 数据分析实验：预测英国的房价。

使用R下载数据文件，因为该数据是不断更新的，写博客当天的时间跨度是91年至17年5月的数据。

library(gdata)
url="http://www.nationwide.co.uk/~/media/MainSite/documents/about/house-price-index/downloads/monthly.xls"
download.file(url,"monthly_jzh.xls",mode='wb')#xls是windows下的binary文件，所以需要加一个mode选项
dt=read.xls("monthly_jzh.xls",perl="C:\\Strawberry\\perl\\bin\\perl.exe",stringsAsFactors=F)##需要指定perl路径
library(zoo)
dt$date=as.yearmon(dt$date)
dt_zoo=read.zoo(dt,index.column=1)
summary(dt_zoo);head(dt_zoo);frequency(dt_zoo)

PS：这里下载的数据文件保存在当前工作目录，因为博主比较渣，看所以不能读取原始数据文件，这里在运行read.xls之前

是对下载下来的文件进行过手动修改的，修改的是前两列的格式，第一列修改为日期格式，第二列设置为数值格式，

顺带把第一列命名为“date”。如果有时间，后期再修改一下吧，难以忍受这种手动干预。

以上代码运行完，最后frequency得到12就对了，这实际上是一个英国的月度房价数据，现在是老外的数据拿来做研究

，以后爬取国内的数据来分析试试（比如某家），毕竟房价是个比较热门的话题。

然后我们就可以开始建模了。

2.1.1 第一步：模型识别和参数估计

R语言forecast包中，一个auto.arima函数就可以实现模型识别和参数估计。

library(forecast)
rtn=diff(dt_zoo)/lag(dt_zoo,k=-1)*100 #计算连续复利收益率
arima_fit=auto.arima(rtn,stationary=TRUE,seasonal = FALSE,ic="aic");arima_fit

结果如下：

Series: rtn 
ARIMA(2,0,0)            with non-zero mean 

Coefficients:
         ar1     ar2    mean
      0.2361  0.3378  0.4381
s.e.  0.0528  0.0529  0.1300

sigma^2 estimated as 0.9919:  log likelihood=-445.79
AIC=899.57   AICc=899.7   BIC=914.6

这里的程序运行结果显示，最终选定AR(2)模型。因为这个数据与教材中的数据选取的时间段不一样，所以结果与教

材的不一样，但是通过观察发现：1、都是选择了AR(2)模型，且系数相差不大，说明这种自回归是非常显著的，并且

这种线性关系比较稳定。关于函数的详细参数及说明，以后另外写（抄袭）文章来补充学习一下。

2.1.2 第二步：模型诊断

模型诊断就是要判断模型建立以后，验证集残差是否符合白噪声的特点，可以通过画图和统计检验进行验证。

tsdiag(arima_fit)

以上代码的输出结果如下图所示：

图一的standardized residuals是残差的散点图（柱状化处理），主观判断波动无明显规律性；图二是自相关图（ACF）

，我也不知道为何有小数。虽然显示一阶和二阶还是有相关性的，但是基本可见相关性不大了；最后的图三通过

Ljung-Box统计量检验自相关性，P至都很大，远在0.8以上，可见不能拒绝原假设，即自相关系数为0的假设。

这样模型就建立完毕了，最后我们的模型就定为一个带截距项（mean）的二阶自回归模型，系数如上。为了可视化

一下我们模型拟合的效果，有一个感性的认识 ，可以通过画一个时间序列曲线，看看一致性如何：

plot(arima_fit$x, lty = 1, main = "预测房价与真实房价对比图", ylab = "Return in percent", xlab = "Date")
lines(fitted(arima_fit), lty = 2,lwd = 2, col = "red")

输出结果为下图：

我看了这个之后觉得，这个预测的真的很不好！！！勉强用这个数据模型预测接下来几个月的房价吧，过些日子来

看看预测能力如何：

predict(arima_fit, n.ahead=3)

预测值为：

$pred
     Jun        Jul        Aug
2017 0.3653420  0.4375064  0.4133708

3、协整

本章的总体行文内容是：协整定义为两个或几个不平稳时间序列，其本身各自均不平稳，但是其存在某种线性

组合，这种线性组合是平稳的。这种组合的线性平稳，我们称其为是一种长期平稳关系（long-termrelationships）。

本节对比了不考虑长期平稳关系的情况和考虑协整关系的模型，查看后者的模型拟合效果会好一些。下面通过一个

例子演示这种区别。

在衍生品市场中，经常会用到交叉（cross）对冲，就是用市场上存在的衍生品品种去对冲我的商品价格风险，

衍生品的基础资产可能和我的商品不相同，但是如果相近，也是可以通过配置一定的比率进行对冲。这种相近的证明

和比率的确定，就可以通过协整模型求解。本例子使用的是喷气式飞机煤油（jet fuel）的对冲，使用的衍生品是加热

油（heating oil）远期产品。用到的主要R包是urca。

PS:鉴于找到了原书的代码和数据，以后会直接粘贴代码，自己加上必要的备注了。

library("zoo")
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)#首先使用不带截距项的回归求解对冲比.
#以上代码使用diff函数是明显的，因为风险敞口是价格的波动值，diff计算价格增加值。
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。

下面我们使用协整关系探索长期平稳关系，代码如下：

3.1 单位根检验

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") # 使用augmented Dickey-Fuller检验，检验商品价格的平稳性
summary(jf_adf) # 原假设是存在单位根，即平稳，p-value较小，拒绝原假设。即不平稳。

ho_adf <- ur.df(prices$HeatingOil, type = "drift") # 对衍生品价格序列进行单位根检验
summary(ho_adf) # 检验同样发现不平稳（P值比较大）

3.2 两步法ECM模型

## 两步法的第一步是对原始数据先做一个线性回归得到残差，检验残差是否平稳，也即线性组合是否平稳
mod_static <- summary(lm(prices$JetFuel ~ prices$HeatingOil)) # 原始数据(不是价格变动值)拟合结果
error <- residuals(mod_static) # 确定估计的残差
error_cadf <- ur.df(error, type = "none") # 对残差检验平稳性,注意参数“type”
summary(error_cadf) #查看残差序列是否是平稳序列

## 两步法的第二步对上述得到的残差，加入到线性模型中，也即建立ECM模型
djf <- diff(prices$JetFuel) # 依旧研究的是价格变动
dho <- diff(prices$HeatingOil)  #同上
error_lag <- lag(error, k = -1) #需要加入模型的残差项
mod_ecm <- lm(djf ~ dho + error_lag) # 加入了残差项
summary(mod_ecm) #输出ECM模型的估计结果

输出结果如下：

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

Residuals:
     Min       1Q   Median       3Q      Max 
-0.19332 -0.03391 -0.00101  0.02128  0.44942 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.001511   0.005013   0.301    0.763    
dho          0.899489   0.032547  27.637   <2e-16 ***
error_lag   -0.655301   0.066303  -9.883   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06892 on 187 degrees of freedom
Multiple R-squared:  0.8189,	Adjusted R-squared:  0.817 
F-statistic: 422.9 on 2 and 187 DF,  p-value: < 2.2e-16

可见，ECM模型求得的最佳对冲比率为0.90020，对应的残差标准差为0.06875。比之前稍微效果好一些。

4、波动率模型（异方差）

在ARIMA模型中，我们是假设序列是平稳的，这样每时每刻的波动方差都是相等的。但是现实中，这种假设往往

不存在，所以就产生了波动率模型。波动率模型的结果经常用在VaR模型。这里的VaR模型并不是统计学中的VAR模型

而是Value at Risk的简写。在建立波动率模型的时候，我们首先检验是否存在波动率的不平稳性，然后在不平稳的情况

下建立ARCH或者GARCH模型，在应用于VaR模型。

library("zoo")
intc <- read.zoo("intc.csv", header = TRUE, sep = ",", format = "%Y-%m", FUN = as.yearmon) #读取英特尔的股票收益数据
plot(intc, main = "Monthly returns of Intel Corporation", xlab = "Date", ylab = "Return in percent") # 收益序列图
### 检验是否存在ARCH effect，两种方法：LB检验和LM检验
# LB检验是直接检查表征方差的平方项的自相关性，以下代码进行 Ljung-Box检验
Box.test(coredata(intc^2), type = "Ljung-Box", lag = 12) # 原假设是不存在自相关性，即不存在ARCH effect，P值接近于0，拒绝原假设。
# LM检验使用FinTS的ArchTest函数进行，代码如下
install.packages("FinTS")
library("FinTS")
ArchTest(coredata(intc)) # 结论相同，拒绝原假设

以上的L-B检验和L-M检验都说明，英特尔股票的收益率存在着波动性。以下建立GARCH模型：

library("rugarch")
intc_garch11_spec <- ugarchspec(variance.model = list(garchOrder = c(1, 1)), mean.model = list(armaOrder = c(0, 0))) # 模型识别
intc_garch11_fit <- ugarchfit(spec = intc_garch11_spec, data = intc) # 模型参数估计
intc_garch11_fit # 查看模型结果

这里需要说明，之所以选择GARCH（1,1）模型是因为该模型比较常用。另外均值模型为白噪声，所以为ARMA(0,0)

这里的模型识别，实际上是人为主观设定。设定完以后输入数据跑极大似然估计模型。

模型构建完毕后需要进行回测，以验证模型。这里构建VaR模型预测并与实际收益率进行比较，如果收益率比

VaR损失更大的概率应该在设定的比例之内，这里指的是1%。下面是要用到的回测函数ugarchroll的一些参数设置：

• 使用的收益率数据存储在zoo对象intc中。
• 回测的起始期（n.start）应该是序列开始（就是1983年1月）时的前120个月。
• 每月都需要重新估计模型（refit.every = 1）。
• 我们使用移动（moving）窗口来估计。
• 我们使用一个混合（hybrid）的解决方法。
• 我们希望计算VAR在99%的尾部置信水平上（VaR.alpha = 0.01）的临界值（calculate.VaR = TRUE）。
• 我们希望保留估计的系数（keep.coef = TRUE）。
• 摘自：http://www.epubit.com.cn/book/onlinechapter/52471；著　[匈牙利] Gergely Daróczi 等；译　高　蓉　 李　茂

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(intc_garch11_roll, type = "VaR", VaR.alpha = 0.01, conf.level = 0.99) # 回测结果

回测的report结果与教材的不一样，我也很挠头啊，回头补充学习VaR和GARCH，做一个chapter 1的辅助材料，写完

以后补充链接。

以下使用图像的形式展现模拟效果，对比VaR和真实收益率的走势图：

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

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

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

预测：

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

根据输出结果得到的结论是：

波动率（sigma）的一步预期是0.1168。因为我们假定了正态分布，置信水平为99%的VaR可以使用标准正态分

布的99%分位数（输入qnorm(0.99)）来计算。因此对下一个周期，一个月的99%VaR就是

qnorm(0.99)*0.1168 =0.2717。

结果，月收益率高于−27%的概率是99%。

好，现在第一章基本抄写结束了，可以看得出，不明白VaR模型和GARCH模型，学习这一张的后半部分还是云

里雾里的。接下去会补充VaR和GARCH的理论，可能一篇博客完成，如果理论内容较多，就分成两篇或若干篇。