R Notebook
\[ ARIMA(p,d,q)(P, D, Q)_{m} \] (p,d,q)为非季节部分,(P, D, Q)_{m}为季节部分，m = 每年的观测数量
比如，一个\(ARIMA(1,1,1)(1,1,1)_4\)模型(不含常数)对应的是季节数据(m=4),可以表示为 \[ (1 - \phi_{1}B)~(1 - \Phi_{1}B^{4}) (1 - B) (1 - B^{4})y_{t} = (1 + \theta_{1}B)~ (1 + \Theta_{1}B^{4})\varepsilon_{t}. \] 即季节性项与非季节性项相乘以后得到。
自相关图/偏自相关图观察季节性特征
AR 或者 MA 模型中的季节性特征其实可以从偏自相关图和自相关图中的季节性延迟中被观察到。
例如 ARIMA(0,0,0)(0,0,1)_12模型会显示出: a) 自相关图中延迟期数为 12 上有一个明显突起，但是没有其他明显的突起； b) 偏自相关图中季节性延迟（即当延迟为12,24,36,…）上出现指数性的衰减.
一个ARIMA(0,0,0)(1,0,0)_12： a) 自相关图中季节性延迟上出现指数型衰减； b) 偏自相关图中延迟为 12 上的一个明显突起。
实例：欧洲季度零售贸易额
library(fpp2)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── fpp2 2.4 ──## ✓ ggplot2 3.2.1 ✓ fma 2.4
## ✓ forecast 8.13 ✓ expsmooth 2.3## #library(feasts)
#library(fabletools)
#eu_retail <- as_tsibble(fpp2::euretail)
#eu_retail %>% autoplot(value) + ylab("Retail index") + xlab("Year")
autoplot(euretail) + ylab("零售指数") + xlab("年份") +
theme(text = element_text(family = "STHeiti"))+
theme(plot.title = element_text(hjust = 0.5))
### 季节性差分
#数据显然是非平稳的，并且拥有一些季节性的特征，因此我们先对进行季节性差分，季节性差分数据如图8.17所示。
euretail %>% diff(lag=4) %>% ggtsdisplay()
#由于非平稳性仍然存在，我们继续进行一阶差分，
euretail %>% diff(lag=4) %>% diff() %>% ggtsdisplay()
经过了差分变换，通过自相关图和偏自相关图来找到合适的 ARIMA 模型。 ACF图中在延迟为1时的突起显示出数据包含非季节性 MA(1) 的成分，而延迟为4的突起则显示出季节性 MA(1)的成分。 因此，ARIMA(0,1,1)(0,1,1)_4 模型开始，这表明数据经过了一次季节性差分和一次一阶差分，同时包含季节性和非季节性的 MA(1) 成分。
ndiffs确定一阶季节性差分次数
usmelec_stl <- stl(usmelec,t.window=13, s.window="periodic", robust=TRUE)
#head(usmelec)
usmelec_stl_s <- usmelec_stl$time.series[,"seasonal"]
usmelec_stl_t <- usmelec_stl$time.series[,"trend"]
usmelec_stl_r <- usmelec_stl$time.series[,"remainder"]
F_T <- max(0,1-var(usmelec_stl_r)/var(usmelec_stl_t+usmelec_stl_r))
#F_T
F_S <- max(0,1-var(usmelec_stl_r)/var(usmelec_stl_s+usmelec_stl_r))
F_S## [1] 0.7864134如果F_S<0.64，不需要进行季节性差分，否则需要进行一次季节性差分。
usmelec %>% log() %>% nsdiffs()## [1] 1usmelec %>% log() %>% diff(lag=12) %>% ndiffs()## [1] 1usmelec %>% log() %>% diff(lag=12) %>% diff() %>% ndiffs()## [1] 0由于nsdiffs()函数返回1，说明需要进行一次季节性差分，我们对季节差分后的数据运行ndiffs()函数。这些函数的运行结果说明我们应该进行一次季节性差分和一次一步差分。
确定模型
euretail %>%
Arima(order=c(0,1,1),seasonal=c(0,1,1)) %>%
residuals() %>%
ggtsdisplay()
自相关图和偏自相关图都在延迟为2的地方出现了明显的突起，在延迟为3的地方出现了较为明显的突起，这反映出有一些额外的非季节性项需要被补充进入模型。
fit1 <- euretail %>% Arima(order=c(0,1,2),seasonal=c(0,1,1))
fit2 <- euretail %>% Arima(order=c(0,1,3),seasonal=c(0,1,1))
rbind(fit1$aicc,fit2$aicc)## [,1]
## [1,] 74.26992
## [2,] 68.39031ARIMA(0,1,2)(0,1,1)_4 模型的AICc值 为 74.26，而 ARIMA(0,1,3)(0,1,1)_4 模型的AICc值为 68.39。 因此,我们最后选择了 ARIMA(0,1,3)(0,1,1)_4模型。
预测
fit3 <- Arima(euretail,order=c(0,1,3),seasonal=c(0,1,1))
checkresiduals(fit3)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,3)(0,1,1)[4]
## Q* = 0.51128, df = 4, p-value = 0.9724
##
## Model df: 4. Total lags used: 8所有残差都在显著性临界值内，因此它们类似于白噪声。Ljung-Box 检验也显示残差之间不存在自相关性。
fit3 %>% forecast(h=12) %>% autoplot()
### auto.arima()自动确定模型 我们可以用auto.arima()来完成这里大部分的工作，它的结果如下所示:
fit.arima <- auto.arima(euretail)
fit.arima## Series: euretail
## ARIMA(0,1,3)(0,1,1)[4]
##
## Coefficients:
## ma1 ma2 ma3 sma1
## 0.2630 0.3694 0.4200 -0.6636
## s.e. 0.1237 0.1255 0.1294 0.1545
##
## sigma^2 estimated as 0.156: log likelihood=-28.63
## AIC=67.26 AICc=68.39 BIC=77.65实例：澳大利亚皮质类固醇药物销量
平稳定方差
# 澳大利亚皮质类固醇药物销量（每月每百万剂）
#图中方差随着销量水平增加而略微增加，因此我们运用log来平稳定方差。
library(fpp2)
head(h02)## Jul Aug Sep Oct Nov Dec
## 1991 0.429795 0.400906 0.432159 0.492543 0.502369 0.602652lh02 <- log(h02)
cbind("H02 sales (million scripts)" = h02,
"Log H02 sales"=lh02) %>%
autoplot(facets=TRUE) + xlab("年份") + ylab("")
### 季节性差分 数据拥有很强的季节性，并且显然非平稳，因此我们使用季节性差分 不清楚是否应该再做一次差分，我们决定不继续差分
lh02 %>% diff(lag=12) %>%
ggtsdisplay(xlab="年份",main="季节性差分后的H02药物销量")
在季节性差分后的数据的图中，我们发现偏自相关图在延迟为 12 和 24 的地方存在突起，然而在自相关图中并不存在季节性延迟的突起。 这有可能暗示出模型应含有季节性 AR(2) 的部分。在非季节性的延迟中，偏自相关图有三个明显突起，显示出可能存在的 AR(3) 部分。自相关图并没有任何简单模型的特征，因此我们不作考虑。
确定了数据对应的一个可能的模型 ARIMA(3,0,0)(2,1,0)_12 我们用这个模型以及它的一些衍生模型进行拟合，并且计算出了各自的AICc值
注意： 当我们使用 AICc 值来比较模型时，所有的模型应该是\(\color{red}{同阶差分}\)的，这一点很重要。
确定模型
fit1 <- lh02 %>% Arima(order=c(3,0,1),seasonal=c(0,1,2))
fit2 <- lh02 %>% Arima(order=c(3,0,1),seasonal=c(1,1,1))
fit3 <- lh02 %>% Arima(order=c(3,0,1),seasonal=c(0,1,1))
fit4 <- lh02 %>% Arima(order=c(3,0,1),seasonal=c(2,1,0))
fit5 <- lh02 %>% Arima(order=c(3,0,0),seasonal=c(2,1,0))
fit6 <- lh02 %>% Arima(order=c(3,0,2),seasonal=c(2,1,0))
fit7 <- lh02 %>% Arima(order=c(3,0,1),seasonal=c(1,1,0))
rbind(fit1$aicc,fit2$aicc,fit3$aicc,fit4$aicc,fit5$aicc,fit6$aicc,fit7$aicc)## [,1]
## [1,] -485.4754
## [2,] -484.2511
## [3,] -483.6674
## [4,] -476.3086
## [5,] -475.1230
## [6,] -474.8826
## [7,] -463.3951在这些模型中，ARIMA(3,0,1)(0,1,2)_12模型拥有最小的AICc值，因而是最优模型
Box-Cox transformation parameter. If lambda=“auto”, then a transformation is automatically selected using BoxCox.lambda. The transformation is ignored if NULL. Otherwise, data transformed before model is estimated \[ \begin{align*} &\text{Box-Cox transformation}\\ & y(\lambda) =\frac{y^\lambda - 1}{\lambda} , \lambda != 0 \\ & y(\lambda) =log(y) , \lambda = 0 \\ \end{align*} \]
(fit <- Arima(h02,order=c(3,0,1),seasonal=c(0,1,2),lambda=0))## Series: h02
## ARIMA(3,0,1)(0,1,2)[12]
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ar3 ma1 sma1 sma2
## -0.1603 0.5481 0.5678 0.3827 -0.5222 -0.1768
## s.e. 0.1636 0.0878 0.0942 0.1895 0.0861 0.0872
##
## sigma^2 estimated as 0.004278: log likelihood=250.04
## AIC=-486.08 AICc=-485.48 BIC=-463.28checkresiduals(fit,lag=36)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(3,0,1)(0,1,2)[12]
## Q* = 50.712, df = 30, p-value = 0.01045
##
## Model df: 6. Total lags used: 36auto.arima()自动确定模型
fit_auto <- auto.arima(h02,stepwise = TRUE,seasonal = TRUE,ic='aicc',lambda=0)
fit_auto## Series: h02
## ARIMA(2,1,1)(0,1,2)[12]
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 ar2 ma1 sma1 sma2
## -1.1358 -0.5753 0.3683 -0.5318 -0.1817
## s.e. 0.1608 0.0965 0.1884 0.0838 0.0881
##
## sigma^2 estimated as 0.004278: log likelihood=248.25
## AIC=-484.51 AICc=-484.05 BIC=-465接下来我们将会使用自动的 ARIMA 算法，在默认设置下，auto.arima()函数返回了 ARIMA(2,1,1)(0,1,2)_12模型。然而，这个模型同样没有通过 Ljung-Box 检验。有时候找到一个通过所有检验的模型是不可能的。
模型RMSE]
包含最近两年数据的测试集对之前拟合出的模型进行比较。 这些模型中没有一个能够通过所有的残差检验，在实际应用中，我们仍然会使用我们找到的最优模型，即使它不能通过所有的检验。
train <- window(h02, start=c(1991,7),end=c(2006,6)) test <- window(h02, start=c(2006,7),end=c(2008,6))
a1 <- fit1 %>% forecast(h=12) %>% accuracy()
a2 <- fit2 %>% forecast(h=12) %>% accuracy()
a3 <- fit3 %>% forecast(h=12) %>% accuracy()
a4 <- fit4 %>% forecast(h=12) %>% accuracy()
a5 <- fit5 %>% forecast(h=12) %>% accuracy()
a6 <- fit6 %>% forecast(h=12) %>% accuracy()
a7 <- fit7 %>% forecast(h=12) %>% accuracy()
a8 <- fit_auto %>% forecast(h=12) %>% accuracy()
rbind(a1[2],a2[2],a3[2],a4[2],a5[2],a6[2],a7[2],a8[2])## [,1]
## [1,] 0.06245679
## [2,] 0.06269911
## [3,] 0.06321256
## [4,] 0.06445510
## [5,] 0.06509306
## [6,] 0.06441930
## [7,] 0.06751302
## [8,] 0.05015710ARIMA(2,1,1)(0,1,2) _12模型的预测（拥有最小的均方根误差,并且在只进行季节性差分的模型中拥有最优的AICc值）
h02 %>%
Arima(order=c(2,1,1),seasonal=c(0,1,2),lambda=0) %>%
forecast() %>%
autoplot() +
ylab("H02销量(百万剂)") + xlab("年份")+
theme(text = element_text(family = "STHeiti"))+
theme(plot.title = element_text(hjust = 0.5))
案例（知乎）:
https://zhuanlan.zhihu.com/p/88442615
ARIMA和ETS都是时间序列分析中的有效工具，在实践中经常将两者进行比较。
library(fpp2)
head(qcement)## Qtr1 Qtr2 Qtr3 Qtr4
## 1956 0.465 0.532 0.561 0.570
## 1957 0.529 0.604#cement <- window(qcement, start=1988)
#train <- window(cement, end=c(2007,4))
train <- window(qcement, start=c(1988,1),end=c(2007,4))
(fit.arima <- auto.arima(train))## Series: train
## ARIMA(1,0,1)(2,1,1)[4] with drift
##
## Coefficients:
## ar1 ma1 sar1 sar2 sma1 drift
## 0.8886 -0.2366 0.081 -0.2345 -0.8979 0.0105
## s.e. 0.0842 0.1334 0.157 0.1392 0.1780 0.0029
##
## sigma^2 estimated as 0.01146: log likelihood=61.47
## AIC=-108.95 AICc=-107.3 BIC=-92.63checkresiduals(fit.arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1)(2,1,1)[4] with drift
## Q* = 3.3058, df = 3, p-value = 0.3468
##
## Model df: 6. Total lags used: 9# 进行预测
a1 <- fit.arima %>% forecast(h = 4*(2013-2007)+1) %>% accuracy(qcement)
a1[,c("RMSE","MAE","MAPE","MASE")]## RMSE MAE MAPE MASE
## Training set 0.1001195 0.07988903 4.372443 0.5458078
## Test set 0.1996098 0.16882205 7.719241 1.1534049