1 目的
该篇文章主要展示针对时序进行ARIMA乘法模型建模，并根据实际情况进行改进。案例数据同 时间序列分析实战（三）：时序因素分解法：某欧洲小镇1963年1月至1976年12月每月旅馆入住的房间数构成时间序列x t x_txt，我们在 时间序列分析实战（五）：ARIMA加法（疏系数）模型建模利用此数据构建了ARIMA加法模型。现考虑演示根据该数据建立ARIMA乘法模型。
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乘法模型建立
运行程序：
par(mfrow=c(2,1)) acf(y,lag=36) title(sub="图2 入住房间数差分后序列自相关图(ACF)") pacf(y,lag=36) title(sub="图3 入住房间数差分后序列偏自相关图(PACF)")
运行结果：
图2 入住房间数差分后序列自相关图(ACF)）
图3 入住房间数差分后序列偏自相关图(PACF)
为考虑季节相关特征，这时考察延迟12阶、24阶等以周期长度为单位的自相关系数和偏自相关系数的特征。图2和图3显示,自相关图表现拖尾性质，偏自相关图表现截尾性质，若把偏自相关系数看作1周期截尾，结合前面的差分信息，经过多次尝试后，尝试拟合ARIMA乘法模型:ARIMA(3,1,2)×(1,1,0)12
6 ARIMA乘法模型拟合
运行程序：
fit4=arima(data1,order = c(3,1,2),seasonal = list(order=c(1,1,0),period=12)) fit4
运行结果：
## ## Call: ## arima(x = data1, order = c(3, 1, 2), seasonal = list(order = c(1, 1, 0), period = 12)) ## ## Coefficients: ## ar1 ar2 ar3 ma1 ma2 sar1 ## 0.8229 0.0030 -0.3107 -1.7905 0.7905 -0.4073 ## s.e. 0.1155 0.1036 0.0851 0.1110 0.1102 0.0816 ## ## sigma^2 estimated as 213.2: log likelihood = -640.63, aic = 1295.25
7 ARIMA乘法模型显著性检验
运行程序：
#模型显著性检验 ts.diag(fit4) title(sub="图4 ARIMA(3,1,2)×(1,1,0)_12模型显著性检验")
运行结果：
图4 ARIMA(3,1,2)×(1,1,0)_12模型显著性检验
此时模型ϕ 2 \phi _2ϕ2估计值在两倍标准差内，考虑建立ARIMA乘法疏系数模型：ARIMA((1,3),1,2)×(1,1,0)12
8 ARIMA乘法疏系数模型
运行程序：
#拟合乘法ARIMA模型 fit5=arima(data1,order = c(3,1,2),seasonal = list(order=c(1,1,0),period=12), transform.pars = F,fixed=c(NA,0,NA,NA,NA,NA)) fit5
运行结果：
## ## Call: ## arima(x = data1, order = c(3, 1, 2), seasonal = list(order = c(1, 1, 0), period = 12), ## transform.pars = F, fixed = c(NA, 0, NA, NA, NA, NA)) ## ## Coefficients: ## ar1 ar2 ar3 ma1 ma2 sar1 ## 0.8263 0 -0.3090 -1.7930 0.793 -0.4072 ## s.e. 0.0941 0 0.0611 0.1078 0.107 0.0807 ## ## sigma^2 estimated as 213.2: log likelihood = -640.63, aic = 1293.25
9 ARIMA乘法疏系数模型显著性检验
运行程序：
#模型显著性检验 ts.diag(fit5) title(sub="图5 ARIMA((1,3),1,2)模型显著性检验")
运行结果：
图5 ARIMA(3,1,2)×(1,1,0)_12模型显著性检验
10 ARIMA乘法疏稀疏模型预测
运行程序：
library(forecast) fore2=forecast::forecast(fit5,h=36) fore2$mean
运行结果：
## Jan Feb Mar Apr May Jun Jul ## 1977 827.6433 768.2262 771.2960 860.6892 838.8292 963.9901 1128.1026 ## 1978 853.5056 785.6367 792.7410 886.3393 869.2913 985.3445 1154.0686 ## 1979 875.9526 811.5455 817.0429 908.9636 889.9786 1009.7484 1176.5900 ## Aug Sep Oct Nov Dec ## 1977 1159.3063 892.6915 897.5451 778.5312 884.5493 ## 1978 1178.2927 915.9224 915.4218 804.8587 914.4704 ## 1979 1203.6458 939.5356 941.2077 827.1998 935.3489
运行程序：
plot(fore2,lty=2,sub='图5 入住房间数序列疏系数的ARIMA乘法模型预测效果图') lines(fore2$fitted,col=4)
运行结果：
图6 入住房间数序列疏系数的ARIMA乘法模型预测效果图
根据结果显示，疏稀疏的ARIMA乘法模型有较好拟合效果。
