ARIMA模型
自回归移动平均模型（ARIMA）包含一个确定（explicit）的统计模型用于处理时间序列的不规则部分，它允许不规则部分可以自相关。
数据准备
> library(zoo)
> library(xts)
> library(TTR)
> kings <- scan("http://robjhyndman.com/tsdldata/misc/kings.dat",skip=3)
Read 42 items
> kingstimeseries <- ts(kings)
> plot.ts(kingstimeseries)
数据差分
如果从一个非平稳的时间序列开始， 那么首先就需要做时间序列差分直到得到一个平稳时间序列。如果必须对时间序列做 d 阶差分才能得到一个平稳序列，那么就使用 ARIMA(p,d,q) 模型，其中 d 是差分的阶数。
> kingstimeseriesdiff1 <- diff(kingstimeseries, differences=1)
> plot.ts(kingstimeseriesdiff1)
从一阶差分的图中可以看出，数据已经是平稳的。如果不平稳就需要继续差分。所以对于英国国王去世年龄时间序列选择 ARIMA(p,1,q) 的模型。
寻找 ARIMA 模型
通常我们需要检查平稳时间序列的自相关图和偏相关图，也就是寻找 ARIMA(p,d,q) 中合适的 p 值和 q 值。
> acf(kingstimeseriesdiff1, lag.max=20)
> acf(kingstimeseriesdiff1, lag.max=20, plot=F)
Autocorrelations of series ‘kingstimeseriesdiff1’, by lag
0 1 2 3 4 5 6 7 8 9 10
1.000 -0.360 -0.162 -0.050 0.227 -0.042 -0.181 0.095 0.064 -0.116 -0.071
11 12 13 14 15 16 17 18 19 20
0.206 -0.017 -0.212 0.130 0.114 -0.009 -0.192 0.072 0.113 -0.093
可以看出在滞后 1 阶的自相关值超出了置信边界，但是其他所有在滞后1-20 阶的自相关值都没有超出置信边界。
> pacf(kingstimeseriesdiff1, lag.max=20)
> pacf(kingstimeseriesdiff1, lag.max=20, plot=F)
Partial autocorrelations of series ‘kingstimeseriesdiff1’, by lag
1 2 3 4 5 6 7 8 9 10 11
-0.360 -0.335 -0.321 0.005 0.025 -0.144 -0.022 -0.007 -0.143 -0.167 0.065
12 13 14 15 16 17 18 19 20
0.034 -0.161 0.036 0.066 0.081 -0.005 -0.027 -0.006 -0.037
显示在滞后 1，2 和 3 阶时的偏自相关系数超出了置信边界，为负值，且在等级上随着滞后阶数的增加而缓慢减少。
> library(forecast)
> kings.arima <- arima(kingstimeseries, order=c(1, 2, 3))
> summary(kings, arima)
Min. 1st Qu. Median Mean 3rd Qu. Max.
13.00 44.00 56.00 55.29 67.75 86.00
模型预测
> kingsarimaforecast <- forecast(kings.arima, h=5)
> kingsarimaforecast
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
43 66.86843 47.05546 86.68139 36.56711 97.16974
44 69.46189 48.05746 90.86632 36.72664 102.19714
45 69.33508 47.11732 91.55285 35.35594 103.31423
46 69.79963 46.62460 92.97466 34.35648 105.24278
47 70.13562 46.03435 94.23689 33.27591 106.99533
> plot(kingsarimaforecast)
模型误差检验
在指数平滑模型下， 观察 ARIMA 模型的预测误差是否服从零均值、方差不变的正态分布，同时也可以观察连续预测误差是否（自）相关。 同时对国王去世年龄使用 ARIMA(1,1,3) 模型后所产生的预测误差做（自）相关图，和 Ljung-Box 检验
acf(kingsarimaforecast$residuals, lag.max=20)
> Box.test(kingsarimaforecast$residuals, lag=20, type="Ljung-Box")
Box-Ljung test
data: kingsarimaforecast$residuals
X-squared = 10.141, df = 20, p-value = 0.9655
ARIMA 模型自动预测
> kings.arima1 <- auto.arima(kingstimeseries)
> summary(kings.arima1)
Series: kingstimeseries
ARIMA(0,1,1)
Coefficients:
ma1
-0.7218
s.e. 0.1208
sigma^2 estimated as 236.2: log likelihood=-170.06
AIC=344.13 AICc=344.44 BIC=347.56
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.9712931 14.99836 11.92162 -10.40664 29.5176 0.7496724 0.05350284
> kings.arima_new <- arima(kingstimeseries, order=c(0, 1, 1))
> summary(kings.arima_new)
Call:
arima(x = kingstimeseries, order = c(0, 1, 1))
Coefficients:
ma1
-0.7218
s.e. 0.1208
sigma^2 estimated as 230.4: log likelihood = -170.06, aic = 344.13
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.9712931 14.99836 11.92162 -10.40664 29.5176 0.7496724 0.05350284
> kingsarimaforecast <- forecast(kings.arima_new, h=5)
> kingsarimaforecast
Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
43 67.75063 48.29647 87.20479 37.99806 97.50319
44 67.75063 47.55748 87.94377 36.86788 98.63338
45 67.75063 46.84460 88.65665 35.77762 99.72363
46 67.75063 46.15524 89.34601 34.72333 100.77792
47 67.75063 45.48722 90.01404 33.70168 101.79958
> plot(kingsarimaforecast)
> acf(kingsarimaforecast$residuals, lag.max=20)
> Box.test(kingsarimaforecast$residuals, lag=20, type="Ljung-Box")
Box-Ljung test
data: kingsarimaforecast$residuals
X-squared = 13.584, df = 20, p-value = 0.8509
指数平滑模型
它对时间序列上面连续的值之间相关性没有要求，指数平滑法可以用于时间序列数据的短期预测。
简单指数平滑法
适用于没有季节性变化且处于恒定水平以及没有明显趋势的时间序列的预测。
获取数据（数据来源为伦敦每年降雨量），通过 ts 函数转换为时间序列.。
> rain <- scan("http://robjhyndman.com/tsdldata/hurst/precip1.dat", skip=1)
Read 100 items
> rainseries <- ts(rain, start=c(1813))
> plot.ts(rainseries)
> rainseriesHW <- HoltWinters(rainseries, beta=F, gamma=F)
> plot(rainseriesHW)
预测未来 5 年的降水量
> library(forecast)
> rainseriesforecasts <- forecast(rainseriesHW, h=5)
> plot(rainseriesforecasts)
蓝线是预测 1913-1920 间的降雨量，深灰色阴影区域为 80% 的预测区间，浅灰色阴影区域为 95% 的预测区间。 forecast 提供了预测误差的统计指标(residuals),来评估预测是否有改进的可能性：如果预测误差是相关的，则很可能是简单指数平滑预测可以被另外一种预测技术优化。
> acf(rainseriesforecasts$residuals[!is.na(rainseriesforecasts$residuals)], lag.max=20, plot=T)
可以发现自相关系数在第 3 期的时候达到了置信界限。为了验证在滞后 1-20 阶时非 0 自相关属性是否显著，可以借助 Box.test() 的 Ljung-Box检验。
> Box.test(rainseriesforecasts$residuals, lag=20, type="Ljung-Box")
Box-Ljung test
data: rainseriesforecasts$residuals
X-squared = 17.401, df = 20, p-value = 0.6268
统计量为 17.4，并且 P 值是 0.626 这样的值不足以拒绝预测误差在 1-20 阶是非零自相关。
霍尔特指数平滑法
霍尔特指数平滑法可以用于非恒定水平，没有季节性可相加模型的时间序列预测。 霍尔特指数平滑法是估计当前时间的水平和斜率。其平滑水平是由两个参数控制，alpha：估计当前点水平，beta：估计当前点趋势部分斜率。两个参数都介于 0-1 之间，当参数越接近 0，大部分近期的观测值的权值将较小。 数据来源是 1866 年到 1911 年每年女士裙子直径，将数据通过 ts 函数转换为时间序列，并画出时序图。
> skirts <- scan("http://robjhyndman.com/tsdldata/roberts/skirts.dat",skip=5)
Read 46 items
> skirtsseries <- ts(skirts, start=c(1866))
> plot.ts(skirtsseries)
> skirtsseriesHW <- HoltWinters(skirtsseries, gamma=F)
> skirtsseriesHW
Holt-Winters exponential smoothing with trend and without seasonal component.
Call:
HoltWinters(x = skirtsseries, gamma = F)
Smoothing parameters:
alpha: 0.8383481
beta : 1
gamma: FALSE
Coefficients:
[,1]
a 529.308585
b 5.690464
> plot(skirtsseriesHW)
相关预测值中 alpha 值为 0.8383，beta预测值为 1.0，这些都是非常高的值，充分显示了无论是水平上还是趋势的斜率上，当前值对时间序列上的最近的观测值的依赖关系比较重，这样的结果也符合我们的预期，因为时间序列的水平和斜率在整个时间段内发生了巨大的变化。 总体来看，预测的效果也还不错（红色为预测值）。
预测未来 5 年的数据值，并画出预测结果。
> skirtsseriesforecasts <- forecast(skirtsseriesHW, h=5)
> plot(skirtsseriesforecasts)
为了检验预测效果，我们同样检验延迟 1-20 阶中的预测误差是否非零自相关，继续采用 Ljung-Box 检验。
> acf(skirtsseriesforecasts$residuals[!is.na(skirtsseriesforecasts$residuals)], lag.max=20, plot=T)
相关图呈现样本内预测误差在滞后 5 阶时超过置信边界，其他都没有超过，我们认为存在一定的偶尔因素。
> Box.test(skirtsseriesforecasts$residuals, lag=20, type="Ljung-Box")
Box-Ljung test
data: skirtsseriesforecasts$residuals
X-squared = 19.731, df = 20, p-value = 0.4749
p =0.4749，意味着置信度只有 53% 这样的值不足以拒绝预测误差在 1-20 阶是非零自相关，则我们接受预测误差在 1-20 阶是非零自相关的。
Holt-Winters 指数平滑法
有增长或者降低趋势并且存在季节性波动的时间序列的预测方法。 Holt-Winters 算法中提供了 alpha、beta 和 gamma 来分别对应当前点的水平、趋势部分和季节部分，参数的去执法范围都是 0-1 之间，并且参数接近 0 时，近期的观测值的影响权重就越小。 数据来源是澳大利亚昆士兰州海滨纪念商品的月度销售日子做为分析对象，将数据通过 ts 函数转换为时间序列，并画出时序图。
> souvenir <- scan("http://robjhyndman.com/tsdldata/data/fancy.dat")
Read 84 items
> souvenirtimeseries <- ts(souvenir, frequency=12, start=c(1987, 1))
> plot.ts(souvenirtimeseries)
可以通过取对数来减少极值带来的影响，消除方差不齐。
> logsouvenirtimeseries <- log(souvenirtimeseries)
> plot.ts(logsouvenirtimeseries)
> souvenirtimeseriesHW <- HoltWinters(logsouvenirtimeseries)
> plot(souvenirtimeseriesHW)
通过 forecast 包来预测未来 12 个月的销售数据，并画出预测结果
> souvenirtimeseriesforecasts <- forecast(souvenirtimeseriesHW, h=12)
> plot(souvenirtimeseriesforecasts)
模型非常成功得预测了季节峰值，峰值大约发生在每年的 12 月份。 还可以通过画相关图和进行 Ljung-Box 检验来检查样本内预测误差在延迟 1-20 阶时否是非零自相关的，并以此确定预测模型是否可以再被优化。
> acf(souvenirtimeseriesforecasts$residuals[!is.na(souvenirtimeseriesforecasts$residuals)], lag.max=20, plot=T)
相关图显示出在滞后 1-20 阶中样本自相关值都没有超出显著（置信）边界。
> Box.test(souvenirtimeseriesforecasts$residuals, lag=20, type="Ljung-Box")
Box-Ljung test
data: souvenirtimeseriesforecasts$residuals
X-squared = 17.53, df = 20, p-value = 0.6183
Ljung-Box 检验的 p 值为 0.6183，所以我们推断在滞后 1-20 阶中没有明显证据说明预测误差是非零自相关的。
