作者：黄天元，复旦大学博士在读，热爱数据科学与开源工具（R），致力于利用数据科学迅速积累行业经验优势和科学知识发现，涉猎内容包括但不限于信息计量、机器学习、数据可视化、应用统计建模、知识图谱等，著有《R语言数据高效处理指南》（《R语言数据高效处理指南》(黄天元)【摘要 书评 试读】- 京东图书，《R语言数据高效处理指南》(黄天元)【简介_书评_在线阅读】 - 当当图书）。知乎专栏：R语言数据挖掘。邮箱：huang.tian-yuan@qq.com.欢迎合作交流。
ARIMA中加入新的信息进行回归
在上一部分（HopeR：R语言时间序列分析（十二）：ARIMA）中讨论了ARIMA的使用，但是之前提到的所有时间序列预测基本都只用了自己来预测自己，没有额外的任何信息。因此我们所有的假设就是事物的发展具有一定的趋势，具有一定的周期性，新近的状态对未来预测的效果更佳，随着时间的变化周期性和趋势性也有变化等等。但其实，一个事物的发展，往往会受到其他因素的影响。比如个人的消费当然具有一定的时间规律，比如每月几号要交租等等，但是根本上来讲它也与个人的收入有关系（没有钱挣怎么有钱花）。ARIMA允许把一个额外的一个或多个变量纳入其中，然后作为额外的补充信息，进行回归，演示如下：
library(fpp2)
(fit <- auto.arima(uschange[,"Consumption"],
xreg=uschange[,"Income"]))
#> Series: uschange[, "Consumption"]
#> Regression with ARIMA(1,0,2) errors
#>
#> Coefficients:
#> ar1 ma1 ma2 intercept xreg
#> 0.692 -0.576 0.198 0.599 0.203
#> s.e. 0.116 0.130 0.076 0.088 0.046
#>
#> sigma^2 estimated as 0.322: log likelihood=-156.9
#> AIC=325.9 AICc=326.4 BIC=345.3这里，xreg参数可以放入额外的解释变量，可以看到其回归系数xreg为0.203。公式的解读见9.2 Regression with ARIMA errors in R | Forecasting: Principles and Practice，然后让我们对其残差进行诊断：
checkresiduals(fit)
Ljung-Box test
data: Residuals from Regression with ARIMA(1,0,2) errors
Q* = 5.8916, df = 3, p-value = 0.117
Model df: 5. Total lags used: 8
利用上述模型进行预测
需要注意的是，如果要用这个模型进行预测，那么就必须要有未来这个变量的信息。举个例子就是，你想知道你未来支出是多少，你也要对你未来挣多少钱有一个预估。具体实现如下：
fcast <- forecast(fit, xreg=rep(mean(uschange[,2]),8)) #假设未来8个时段都是均值
autoplot(fcast) + xlab("Year") +
ylab("Percentage change")
多变量回归
上面是单变量回归，多元回归是类似的，把变量一个一个键入即可，如下：
xreg <- cbind(MaxTemp = elecdaily[, "Temperature"],
MaxTempSq = elecdaily[, "Temperature"]^2,
Workday = elecdaily[, "WorkDay"])
fit <- auto.arima(elecdaily[, "Demand"], xreg = xreg)
fit
Series: elecdaily[, "Demand"]
Regression with ARIMA(2,1,2)(2,0,0)[7] errors
Coefficients:
ar1 ar2 ma1 ma2 sar1 sar2 drift MaxTemp MaxTempSq Workday
-0.0622 0.6731 -0.0234 -0.9301 0.2012 0.4021 -0.0191 -7.4996 0.1789 30.5695
s.e. 0.0714 0.0667 0.0413 0.0390 0.0533 0.0567 0.1091 0.4409 0.0084 1.2891
sigma^2 estimated as 43.72: log likelihood=-1200.7
AIC=2423.4 AICc=2424.15 BIC=2466.27加入傅里叶项进行回归
如果能够确定在短期内有一定的固定周期性，但是其周期性又是未知的，那么加入傅里叶项将会是一个很好的选择。
cafe04 <- window(auscafe, start=2004)
plots <- list()
for (i in seq(6)) {
fit <- auto.arima(cafe04, xreg = fourier(cafe04, K = i),
seasonal = FALSE, lambda = 0)
plots[[i]] <- autoplot(forecast(fit,
xreg=fourier(cafe04, K=i, h=24))) +
xlab(paste("K=",i," AICC=",round(fit[["aicc"]],2))) +
ylab("") + ylim(1.5,4.7)
}
gridExtra::grid.arrange(
plots[[1]],plots[[2]],plots[[3]],
plots[[4]],plots[[5]],plots[[6]], nrow=3)
请选择AICC值最小的模型。
滞后预测因子
如果你在高考期间，参加了一次又一次的模拟考试，怎么才能预测你下一次模拟考试的成绩？
有人说，用上一次的成绩作为重要参考。
有人说，用上两次的成绩作为重要参考。
有人说，用上三次的成绩作为重要参考。
有人说，用上四次的成绩作为重要参考。
好的不要争了，我们都试试。将成绩的1-3阶滞后项作为预测因子，放入回归即可。
# Lagged predictors. Test 0, 1, 2 or 3 lags.
Advert <- cbind(
AdLag0 = insurance[,"TV.advert"],
AdLag1 = stats::lag(insurance[,"TV.advert"],-1),
AdLag2 = stats::lag(insurance[,"TV.advert"],-2),
AdLag3 = stats::lag(insurance[,"TV.advert"],-3)) %>%
head(NROW(insurance))
# Restrict data so models use same fitting period
fit1 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1],
stationary=TRUE)
fit2 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:2],
stationary=TRUE)
fit3 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:3],
stationary=TRUE)
fit4 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:4],
stationary=TRUE)
c(fit1[["aicc"]],fit2[["aicc"]],fit3[["aicc"]],fit4[["aicc"]])
#> [1] 68.500 60.024 62.833 65.457fit2最棒，因此使用自身和1阶滞后项一起最能够模拟出时间序列的变化模式。
参考资料：
https://otexts.com/fpp2/dynamic.html原文作者：；如有侵权，请联系管理员删除。
