嗯,本来是要学习《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的理论,可能一篇博客完成,如果理论内容较多,就分成两篇或若干篇。