学习资料：
http://www.math.pku.edu.cn/teachers/lidf/course/fts/ftsnotes/html/_ftsnotes/fts-ltscases-gas.html
这一章用三个实例来详细讲解如何用R语言和线性时间序列模型分析实际数据， 并展现线性时间序列模型的适用性与局限性。
数据为：
- 1997-01-06到2010-09-27的美国普通汽油价格周数据；
- 1880年1月到2010年8月全球温度异常值的月度数据；
- 美国失业率月度数据，包括首次申领失业救济金人数的序列以及不包括的序列。
这些数据是持续更新的， 也反映了全球或美国经济的重要方面， 其建模问题有足够的代表性。
用时间序列分析或者统计方法建模时， 最常遇到的困难是如何选取一个适当的模型。当数据之间的动态相依性很复杂时， 模型的形式难以确定；当有多个模型都表现很好时， 模型难以选择。
George Box教授关于建模问题有一句名言(Box 1976):
所有的模型都是错误的，但是其中有一些是有用的。
不同的研究目的可能偏向于不同的模型。
时间序列数据建模的一些指导原则：
数据仅是可利用信息的一部分， 专业知识、常识、历史事件等都是需要考虑的可利用信息。
多个模型可能表现相近， 这时并没有一个“正确的”模型， 选择一个就可以。
在预测时，可以结合多个模型来改善预测效果。
建模的过程是从最简单的模型到逐步复杂， 千万不能以为理论上越复杂、理解和掌握的人数越少的模型才是越好的模型。
模型应尽可能选择更简洁的模型， 如果两个模型的表现相近， 一定要选择更简单的一个。这也是避免过度拟合的要求。过度拟合会导致模型的外推预测能力丧失。
1 数据读入与探索性分析
原油价格和汽油价格对美国经济的重要影响：
20实际70年代初期石油危机展示了石油资源的战略意义。
2008年汽油价格高涨影响了居民生活：
交通成本上涨
供热成本上涨
食品，服务价格上涨==> 通货膨胀
降低了消费者在其它消费项目上的可支配收入，可能导致经济萧条
研究美国汽油价格以及原油价格对其影响。数据为美国普通石油价格周数据，从1997-01-06到2010-09-27。时间是每周三, 共717个时间点; 美国原油价格周数据，从1997-01-03到2010-09-24。时间是每周日。共717个时间点。原油价格比汽油价格早3天发布。
读入数据：石油价格和汽油价格分成了两个输入文件。石油价格的输入文件w-petroprice.txt又是一个用单个空格、多个空格、制表符、空格加制表符分隔的文件。这种文件千万不要自己制造出来。目前readr包还不支持这样的文件， 用R原来的read.table()函数来读入：
rm(list = ls())
如果希望用readr包来读入， 最简单的解决办法是用文本编辑器编辑，将制表符都替换成单个空格；下面的R程序用了R语言的文件操作和正则表达式进行替换：
la <- readLines("w-petroprice.txt")
汽油价格对数值的时间序列图：
plot(
图1: 汽油价格对数值序列
原油价格对数值的时间序列图：
plot(
图2: 原油价格对数值序列
将对数汽油价格与对数原油价格同时绘制：
xts.pgs.pus <- xts(cbind(pgs, pus), index(xts.pgs))
图3: 汽油、原油价格对数值序列
为什么要用对数价格而不用原始价格？请看原始价格的图形：
xts.orig.pgs.pus <- xts(cbind(exp(pgs), exp(pus)), index(xts.pgs))
图4: 汽油、原油价格原始序列
从图形可以看出， 原始价格随时间而振幅变大， 两个序列画在同一序列很难看出有同步变化。对数变换是最常见的克服指数增长、将比例关系转换成线性关系的变换。
对数汽油价格对对数原油价格的散点图：
plot(pus, pgs, xlab="ln(Petrol Price)", ylab="ln(Gasoline Price)",
图5: 对数汽油对对数原油价格
可以看出汽油价格与同期原油价格相关性很强，相关系数为：
> cor.test(pus, pgs)
相关系数为0.98。
因为对数价格呈现出缓慢增长随机趋势， 所以使用其差分值（价格的对数变化率）作为建模数据。
汽油价格对数增长率的时间序列图：
plot(
图6: 汽油价格对数增长率序列
汽油价格对数增长率的直方图：
hist(dpgs, xlab="Gasoline Price Log Increase Rate")
图7: 汽油价格对数增长率分布
汽油价格对数增长率的ACF：
acf(dpgs, main="")
图8: 汽油价格对数增长率ACF
明显非白噪声，也不是低阶MA的典型形状。ACF衰减很快，符合线性时间序列要求。
汽油价格对数增长率的PACF：
pacf(dpgs, main="")
图9: 汽油价格对数增长率PACF
PACF在滞后1到5都是显著的，在滞后19处也比较显著。用低阶AR有一定合理性。PACF衰减很快，符合线性时间序列要求。
对数汽油价格序列的单位根检验：
> fUnitRoots::adfTest(pgs, lags=5, type="c")
检验结果是有单位根。但是， 如果扣除一个非随机线性趋势项检验，就没有单位根：
> fUnitRoots::adfTest(pgs, lags=5, type="ct")
下面还是按照对数汽油价格是单位根序列来建模。
2 AR(5)模型
从时间序列图6来看，在2005年有一个高涨：
plot(
图10: 汽油价格对数增长率序列2005年
在2008年有一个大跌：
plot(
图11: 汽油价格对数增长率序列2008年
这些极端值为建模增加了困难。
使用AIC对AR模型定阶：
> ar(dpgs, method="mle")
AIC选择5阶， 这也与PACF的观察结果吻合。
先建立一个含有常数项的AR(5)模型， 再考虑常数项是否需要：
> arima(dpgs, order=c(5,0,0), include.mean=TRUE)
从intercept的估计值与SE的比较来看， 均值是不显著的。（如果估计值落入正负二倍SE范围就不显著）。所以改为一个不含常数项的AR(5)模型：
> resm3 <- arima(dpgs, order=c(5,0,0), include.mean=FALSE); resm3
逐个看AR系数的显著性， 可以发现ar4明显地不显著， ar2也接近于不显著。删除第4个AR系数：
> resm4 <- arima(
修改后AIC值略有降低。这样限制稀疏模型有可能破坏稳定性条件，所以用polyroot()函数求特征根， 看是否都在单位圆外：
> abs(polyroot(c(1, -coefficients(resm4))))
特征根都在单位圆外且距离单位圆较远， 说明线性时间序列模型是合适的。
如果进一步去掉滞后2系数：
resm5 <- arima(
> resm5 <- arima(
AIC值比仅去掉滞后4的模型变大了。所以应选择仅去掉滞后4的系数。
对此模型做诊断图：
tsdiag(resm4, gof=20)
图12: 候选模型1诊断
从诊断图看，标准化残差仍有较大的异常值。ACF已经基本符合白噪声要求。第三幅各个Ljung-Box白噪声检验均不显著， 其中横坐标是检验所用的自相关系数个数，零假设为符合白噪声要求。这些诊断支持了候选模型1。
3 ARMA(1,3)模型
考虑汽油价格序列对数增长率的其它线性时间序列模型。注意到PACF中滞后1最显著，其它显著位置是刚刚超出0.05水平界限。所以考虑用AR(1)；但是，因为ACF和PACF并不在低阶截尾， 所以单用AR或者单用MA可能不够， 尝试用ARMA。
先看AR(1)的表现：
resm6 <- arima(
图13: AR(1)的残差ACF
AR(1)的ACF在滞后3有一个显著非零。尝试ARMA(1,3):
resm7 <- arima(
> resm7 <- arima(
从系数估计值与SE的比较可以看出MA的系数1、系数2都不显著。先去掉MA系数2， 拟合稀疏系数模型：
> resm8 <- arima(
这个模型的AIC比候选模型1略差一点。MA系数1显著性在两可之间， 所以不继续精简。
称此模型为候选模型2（MC2）。
对候选模型2作诊断图：
tsdiag(resm8, gof=20)
标准化残差图仍有大异常值。ACF已经符合白噪声表现；Ljung-Box白噪声检验基本符合白噪声要求。
因为候选模型1和候选模型2是类似的模型， 候选模型1的AIC较优， 所以两者相比选择候选模型1。
4 固定线性趋势模型
考虑用非随机线性趋势对汽油价格对数值序列建模。
tmp.t <- seq(717)
> summary(resm9)
回归结果显著。考察残差序列：
plot(residuals(resm9), type="l", xlab="Time", ylab="Residual")
图14: 非随机线性趋势模型残差序列
序列呈现出一定的缓慢随机水平变化， 提示非平稳。残差的ACF图：
acf(residuals(resm9), main="", lag.max=20)
图15: 非随机线性趋势模型残差的ACF
说明残差序列的ACF是缓慢衰减的，不太适用线性时间序列加非随机线性趋势作为模型。
5 引入石油价格解释变量的模型
考虑用石油价格对数值作为解释变量， 对汽油价格对数值序列建模。使用带有误差自相关的回归模型。
先做简单回归：
resm10 <- lm(dpgs ~ dpus)
> summary(resm10)
可以看出截距项不显著。尽管因为残差可能有自相关使得检验不一定精确， 因为两个序列都是差分序列，原始序列相关系数很大， 所以回归模型没有常数项是正常的。将上述回归模型改为不带截距项的模型：
resm11 <- lm(dpgs ~ -1 + dpus)
> summary(resm11)
查看回归残差的ACF:
acf(resm11$residuals, main="")
图16: 回归模型的残差序列的ACF
查看回归残差的PACF:
pacf(resm11$residuals, main="")
图17: 回归模型的残差序列的PACF
回归残差的ACF和PACF都是快速衰减的， 但都不截尾。从PACF看， 基本在第5阶之后截尾了。
用AIC对残差的AR模型定阶：
ar(resm11$residuals, method="mle")
> ar(resm11$residuals, method="mle")
AIC选择6阶。使用arima()函数估计AR(6):
resm12 <- arima(resm11$residuals, order=c(6,0,0))
> resm12
从估计结果看，均值（intercept）不显著， 滞后6系数不够显著， 滞后4系数不显著。先拟合一个不带常数项的AR(5):
resm13 <- arima(resm11$residuals, order=c(5,0,0), include.mean=FALSE)
> resm13
这个精简的模型的AIC实际上改善了。因为滞后4系数不显著，做一个稀疏系数AR(5):
resm14 <- arima(resm11$residuals, order=c(5,0,0),
> resm14
AIC继续有改善。将外生回归变量加入模型中，直接建立带有自相关误差项的回归模型：
resm15 <- arima(
> resm15
模型的诊断图形：
tsdiag(resm15, gof=20)
图18: 候选模型3的诊断图形
从诊断图看，除了标准化残差还有大异常值以外， 残差的白噪声性质已经基本确认。说明候选模型3是合适的。
将候选模型1与候选模型3进行比较。模型3的模型方差估计降低了23%， AIC也有明显改善。后面可以看到模型3的预测效果也更好。
6 使用滞后石油价格解释变量的模型
AIC变差，所以退回到上一个模型，以AR(9)去掉滞后4、6、7系数的模型为候选模型4（MC4）。注意，这个模型肯定不如模型3， 因为作为自变量的原油价格比汽油价格的发布时间早了10天。
该部分未能完成，rese16 <- lm(dpgs1 ~ -1 + dpus1)未能成功进行检验。
7 样本外预测
这样， 每个序列有715个观测， 对最后的400个做超前一步预测， 最初仅使用315个观测建模。
dpuslag1 <- dpus[-length(dpus)]
[1] 0.02171326 0.01925884 0.02169466
mafe
[1] 0.01537896 0.01285303 0.01554937
结果的超前一步预测根均方误差和平均绝对预测误差列表如下：
tab1 <- data.frame(
| 模型 | RMSFE | MAFE |
|---|---|---|
| M1 | 0.02171 | 0.01538 |
| M2 | 0.01926 | 0.01285 |
| M3 | 0.02169 | 0.01555 |
看一看同期的汽油价格对数增长率的绝对值的分布：
(该部分来自原讲义，summary的读入数据不对)
summary(abs(dpgs1[nmin:nmax]))
## Min. 1st Qu. Median Mean 3rd Qu. Max.## 0.000000 0.006915 0.014015 0.018875 0.025504 0.162002
有一半的对数增长率在0.014以下， 但是最好的M2模型的预报的平均绝对误差也达到0.013。这提示模型预报效果不够理想。
三个模型比较，利用提前3天的石油价格的模型M2预测效果最好， 不利用石油价格的M1和利用10天前石油价格的M3效果相近， 说明利用10天前的石油价格作用不大。
M1模型的预报图，黑色为真实值， 红色为预报值：
xts.p1 <- xts(cbind(dpgs1[nmin:nmax], fcst[nmin:nmax,1]),
图19: 汽油价格对数增长率用M1模型做一步预测
M2模型的预报图，黑色为真实值， 红色为预报值：
xts.p2 <- xts(cbind(dpgs1[nmin:nmax], fcst[nmin:nmax,2]),
图20: 汽油价格对数增长率用M2模型做一步预测
如果从图形看， 可以看出超前一步预测还是有一定效果的。
参考文献
Box, G.E.P. 1976. “Science and Statistics.” Journal of American Statistician Association, 33526–33536.
