对英国房屋价格建模并预测 ---《量化金融R语言初级教程》

时间序列分析

一、实验介绍

1.1 实验内容

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

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

1.2 课程来源

本课程基于 异步社区 的《量化金融 R 语言初级教程》 第1章制作，感谢 异步社区 授权实验楼发布。如需系统的学习本书，请购买《量化金融 R 语言初级教程》。

为了保证可以在实验楼环境中完成本次实验，我们在原书内容基础上补充了一系列的实验指导，比如实验截图，代码注释，帮助您更好得实战。

如果您对于实验有疑惑或者建议可以随时在讨论区中提问，与同学们一起探讨。

1.3 实验知识点

• 时间序列
• ARMA 模型
• 协整关系
• 改进基本的最小方差对冲比

1.4 实验环境

• R version 3.4.1
• Xfce终端

1.5 适合人群

本课程难度适中，属于初级别课程，需要一定基础的 R 和相关的金融背景 。

1.6 代码及数据的获取

**本实验所需的数据，我们已经帮助读者下载下来上传到了实验平台上**

wget http://labfile.oss.aliyuncs.com/courses/882/data.tar.gztar -zxvf data.tar.gz

注：由于实验需要联网环境，所以 非会员 需要在本实验平台上下载相应的 R 安装包，会员 则无需进行这一步操作,直接进入实验内容。

wget http://labfile.oss.aliyuncs.com/courses/882/tarfile.tar.gztar -zxvf tarfile.tar.gz

通过下面的命令依次将 'forecast' 包的相关依赖包进行本地安装，最后再下载安装 forecast 包。

依赖包：'RcppArmadillo'、 'fracdiff'、 'timeDate'、 'lmtest'、 'tseries'、 'nnet'

'urca' 包的依赖包：'quadprog'

sudo R CMD INSTALL /home/shiyanlou/*.tar.gz # 这里的*号表示包的文件名称，已在上图中显示出来了

二、实验内容

进入 R 环境：

sudo R

2.1　使用时间序列数据

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

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

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

>library(zoo)

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

小提示　下载示例代码

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1．模型识别。

2．模型估计。

3．模型诊断检验。

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

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

除了 zoo 包，我们还会使用到 forecast 包的一些方法。如果你还没安装它，那就需要运行下面的命令来安装 forecast 包。（非会员之前本地安装了就可以跳过此步）

> install.packages("forecast")

下载时间有点长，请耐心等待哟~

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

> library(forecast)

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

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

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

> frequency(hp)[1] 12

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

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

2.2.1　模型识别和估计

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

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

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

> modSeries: hp_retARIMA(2,0,0) with non-zero meanCoefficients:　　　　 ar1　　ar2 intercept　　　0.2299 0.3491　　0.4345s.e.　0.0573 0.0575　　0.1519sigma^2 estimated as 1.105: log likelihood=-390.97AIC=789.94　 AICc=790.1　 BIC=804.28

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

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

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

2.2.2　模型诊断检查

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

> tsdiag(mod)

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

图 2-1　时间序列的诊断

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

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

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

注意打开第一张图后请勿关掉，因为第二张图是基于第一张图的

输出显示在图 2-2 中。

图 2-2　英国房屋数据

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

> accuracy(mod)ME　　　RMSE　　 MAE　　　MPE　　MAPE　　MASE0.00120 1.0514　 0.8059　 -Inf　 Inf　　 0.792980241

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

2.2.3　预测

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

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

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

> plot(forecast(mod))

2.3　协整

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

航空燃油的交叉对冲

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

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

（非会员之前本地安装了'urca'包就可以直接载入'urca'包）

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

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

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

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

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

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

> summary(simple_mod)Call:lm(formula = diff(prices$JetFuel) ~ diff(prices$HeatingOil) +　　0)Residuals:　　 Min　　　 1Q　Median　　　3Q　　 Max-0.52503 -0.02968 0.00131 0.03237 0.39602Coefficients:　　　　　　　　　　　 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 ' ' 1Residual standard error: 0.0846 on 189 degrees of freedomMultiple R-squared: 0.7257, Adjusted R-squared: 0.7242F-statistic: 499.9 on 1 and 189 DF, p-value: < 2.2e-16

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

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

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

注意打开第一张图后请勿关掉，因为第二张图是基于第一张图的

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

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

Test regression driftCall: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.38086Coefficients:　　　　　 Estimate Std. Error t value Pr(>|t|)(Intercept) 0.03050　　0.02177　 1.401 0.16283z.lag.1　　-0.01441　　0.01271　-1.134 0.25845z.diff.lag　0.19471　　0.07250　 2.686 0.00789 **---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Residual standard error: 0.159 on 186 degrees of freedomMultiple R-squared: 0.04099, Adjusted R-squared: 0.03067F-statistic: 3.975 on 2 and 186 DF, p-value: 0.0204Value of test-statistic is: -1.1335 0.9865Critical values for test statistics:　　　1pct 5pct 10pcttau2 -3.46 -2.88 -2.57phi1　6.52　4.63　3.81

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

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

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

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

得到的检验统计量是 −8.912，而规模为 200 的样本在 1% 的置信水平上的临界值为 −2.85。所以，我们拒绝非平稳的原假设。因此我们发现了两个协整变量并且可以进行第二步，两个协整变量意味着一个误差修正模型（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.45117Coefficients:　　　　　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 ' ' 1Residual standard error: 0.06875 on 188 degrees of freedomMultiple R-squared: 0.8198, Adjusted R-squared: 0.8179F-statistic: 427.6 on 2 and 188 DF, p-value: < 2.2e-16

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

三、小结

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

• 立即购买《量化金融R语言初级教程》