🍉CSDN小墨&晓末:https://blog.csdn.net/jd1813346972

   个人介绍: 研一｜统计学｜干货分享
         擅长Python、Matlab、R等主流编程软件
         累计十余项国家级比赛奖项，参与研究经费10w、40w级横向

1 目的

  该篇文章主要展示针对时序进行ARIMA加法模型建模，并根据实际情况进行改进。案例数据同 时间序列分析实战（三）：时序因素分解法：某欧洲小镇1963年1月至1976年12月每月旅馆入住的房间数构成时间序列$x_t$。

2 原序列差分处理

  从 时间序列分析实战（三）：时序因素分解法一文中可知，该序列具有趋势和季节效应，进行1阶差分提取趋势效应，12步差分提取季节效应。

  运行程序：

#对原数据进行1阶12步差分
y=diff(diff(data1,12))
plot(y,sub='图1 入住房间数差分后序列时序图')

  运行结果：

图1 入住房间数差分后序列时序图

3 差分后序列平稳性检验

  运行程序：

#差分后序列平稳性检验
library(aTSA)
adf.test(y)

  运行结果：

## Augmented Dickey-Fuller Test 
## alternative: stationary 
##  
## Type 1: no drift no trend 
##      lag    ADF p.value
## [1,]   0 -19.56    0.01
## [2,]   1 -11.01    0.01
## [3,]   2 -10.63    0.01
## [4,]   3  -9.08    0.01
## [5,]   4 -10.57    0.01
## Type 2: with drift no trend 
##      lag    ADF p.value
## [1,]   0 -19.50    0.01
## [2,]   1 -10.98    0.01
## [3,]   2 -10.60    0.01
## [4,]   3  -9.05    0.01
## [5,]   4 -10.53    0.01
## Type 3: with drift and trend 
##      lag    ADF p.value
## [1,]   0 -19.44    0.01
## [2,]   1 -10.94    0.01
## [3,]   2 -10.56    0.01
## [4,]   3  -9.01    0.01
## [5,]   4 -10.49    0.01
## ---- 
## Note: in fact, p.value = 0.01 means p.value <= 0.01

4 差分后序列白噪声检验

  运行程序：

#差分后序列纯随机性检验
for(k in 1:2) print(Box.test(y,lag=6*k,type="Ljung-Box"))

  运行结果：

## 
##  Box-Ljung test
## 
## data:  y
## X-squared = 56.87, df = 6, p-value = 1.941e-10
## 
## 
##  Box-Ljung test
## 
## data:  y
## X-squared = 76.751, df = 12, p-value = 1.713e-11

  通过1阶12步差分后序列时序图(图1)显示差分后的序列没有明显趋势和周期特征了，ADF检验结果显示差分后序列平稳，纯随机性检验结果显示差分后序列为非白噪声序列，适合建模。

5 ARIMA模型建立

  考虑建立ARIMA加法模型：

$$\nabla s{\nabla }^d{x}_t=\frac{\Theta (B)}{\varphi (B)}{ϵ}_t$$

6 ARIMA模型定阶

  运行程序：

par(mfrow=c(1,2))
acf(y)
title(sub="图2 入住房间数差分后序列自相关图(ACF)")
pacf(y)
title(sub="图2 入住房间数差分后序列偏自相关图(PACF)")

  运行结果：

图2 入住房间数差分后序列自相关图(ACF)（左）& 偏自相关图(PACF)（右）

  差分后序列的自相关图和偏自相关图(图2)显示，自相关系数和偏自相关系数都显示出拖尾的性质。经过多次尝试之后，直到ARMA(1,6)，模型才通过了显著性检验，再考虑到之前的差分运算，实际上拟合的是加法模型：$ARIMA(1,1,6)\times (0,1,0{)}_{12}$

7 ARIMA模型拟合

  运行程序：

#拟合加法ARIMA模型
fit2=arima(data1,order = c(1,1,6),seasonal = list(order=c(0,1,0),period=12),transform.pars = F)
fit2

  运行结果：

## 
## Call:
## arima(x = data1, order = c(1, 1, 6), seasonal = list(order = c(0, 1, 0), period = 12), 
##     transform.pars = F)
## 
## Coefficients:
##           ar1      ma1      ma2      ma3      ma4      ma5     ma6
##       -0.3848  -0.5833  -0.5015  -0.2828  -0.0814  -0.0647  0.5148
## s.e.   0.1739   0.1663   0.1787   0.0933   0.1072   0.0893  0.0941
## 
## sigma^2 estimated as 221.7:  log likelihood = -645.52,  aic = 1307.03

8 ARIMA模型显著性检验

  运行程序：

#模型显著性检验
ts.diag(fit2)
title(sub="图3 ARIMA(1,(1,12),6)模型显著性检验")

  运行结果：

图3 ARIMA(1,(1,12),6)模型显著性检验

  此时，参数${\theta }_4$和${\theta }_5$的估计值在两倍标准差外，未通过显著性检验，所以考虑拟合疏系数模型ARIMA(1,(1,2,3,6))，再考虑到之前的差分运算，实际上拟合的是加法模型$ARIMA(1,1,(1,2,3,6))\times (0,1,0{)}_{12}$

9 ARIMA加法疏系数模型

  运行程序：

#拟合加法ARIMA模型
fit3=arima(data1,order = c(1,1,6),seasonal = list(order=c(0,1,0),period=12),transform.pars = F,fixed=c(NA,NA,NA,NA,0,0,NA))
fit3

  运行结果：

## 
## Call:
## arima(x = data1, order = c(1, 1, 6), seasonal = list(order = c(0, 1, 0), period = 12), 
##     transform.pars = F, fixed = c(NA, NA, NA, NA, 0, 0, NA))
## 
## Coefficients:
##           ar1      ma1      ma2      ma3  ma4  ma5     ma6
##       -0.4201  -0.6015  -0.5797  -0.2909    0    0  0.4730
## s.e.   0.1958   0.1918   0.2185   0.0997    0    0  0.0982
## 
## sigma^2 estimated as 211.5:  log likelihood = -646.01,  aic = 1304.02

10 ARIMA加法疏系数模型显著性检验

  运行程序：

#模型显著性检验
ts.diag(fit3)
title(sub="图4 ARIMA(1,(1,12),(1,2,3,6))模型显著性检验")

  运行结果：

图4 ARIMA(1,(1,12),(1,2,3,6))模型显著性检验

  此时模型参数估计值均在2倍标准差外，显著非零，且模型残差为白噪声序列，aic = 1304.02。最终$ARIMA(1,1,(1,2,3,6))\times (0,1,0{)}_{12}$模型为：

$$(1-B)(1-B{)}^{12}=\frac{1i0.60B-0.58B^2-0.29B^3+0.47B^6}{1+0.42B}{ϵ}_t,Var({ϵ}_t)=211.5$$