【R语言】时间序列案例：住宅销量预测的乘法季节模型

数据来自澳大利亚莫纳什大学统计学教授Rob Hyndman创建的时间序列数据库（TSDB），数据库涵盖近800个时间序列数据（经济、水文、气象、农业等）；该数据集提供了从1973年到1995年美国新单户住宅的月销售量数据，共包含265条数据信息（1973年1月到1995年1月）

对原始数据作图，如图1。可以大致判断该数据存在着序列相关性，以及季节性的变动，同时没有缺失值以及明显异常的点。

数据初步分析和预处理

（1）正态性

时间序列模型的建立和预测通常是在假设数据是正态分布的情况下进行的。通过对序列的观察可以发现，各时间段下数据变化范围所显示的方差随时间变化有较大的不同。做出数据QQ图可以看出，住宅月销量不服从正态分布的，于是进行Box-Cox转换。通过计算lambda = 0.667，最终得到变换后的序列数据： Yt=1λ(Xtλ−1)Y_t=\frac{1}{λ}(X_t^λ-1)Yt=λ1(Xtλ−1)
比较变换前后的数据序列QQ图以及W检验结果可以看到，相比于原始数据，Box-Cox变换后的数据正态性有较大的改善。

Shapiro-Wilk normality test

data: houses.bc	data: houses
W = 0.99337	W = 0.99104
p-value = 0.2908	p-value = 0.1053

同时序列方差也由145.24明显降低为10.52，因此，Box-Cox变换后的数据更适合我们进行建模。变换后的序列如下图。

（2）序列相关性

做出变换后的序列与其滞后各阶的散点图，可以看出一阶滞后自相关很强，二阶、三阶较弱，之后更高的阶数自相关几乎为0，而到了十二阶，自相关性又变得很强。

做出自相关函数和偏自相关函数图像：

通过绘制转换序列的ACF和PACF图，注意到序列中存在一个周期，因为每隔12个月便出现显著的相关性，意味着季节成分的周期是12，即d=12。

（3）平稳性

对变换后的序列进行单位根检验得到：
Augmented Dickey-Fuller Test
data: houses.bc
Dickey-Fuller = -3.7804
Lag order = 6
p-value = 0.02045
由于原假设是存在单位根，即序列非平稳，而p值<0.05，因此变换后的序列具有平稳性。

模型建立

模型建立过程按照下面流程图开展。均基于Box-cox变换后的序列。

（1）模型识别

由第三部分自相关函数带有周期性，以及偏自相关函数滞后一阶、二阶、十三阶的值较大，因此可对序列进行季节差分和一阶差分。
a) 仅进行季节差分的结果如下：

进行单位根检验，得到该序列为平稳的：
Augmented Dickey-Fuller Test
data: diff(houses.bc, lag = 12)
Dickey-Fuller = -3.2832
Lag order = 6
p-value = 0.07447
alternative hypothesis: stationary
查看ACF、PCF图：

从ACF图可以看出序列与滞后前1-11项的相关性线性递减，可以判断为AR模型；结合PACF可知，序列具有拖尾特征，一阶显著，12阶显著，36阶、48阶较为显著，季节部分可以判断为AR模型。
b) 仅进行一阶差分的结果如下：

与上述分析步骤相同，通过单位根检验(Dickey-Fuller = -8.8681,p-value<0.01)、ACF/PACF图分析，我们知道：该序列仍存在明显的周期性，但有效地消除了一些趋势，使平稳性得到了极大的改善。

c) 一阶差分和季节差分得到的序列结果如下：

通过一阶差分有效消除了趋势分量影响，通过季节差分消除了周期因素影响，该序列稳定(Dickey-Fuller = -7.084 , p-value = 0.01)，ACF/PACF如图。从图中可以看到在ACF中显著的只有滞后12阶，在PACF中显著的有滞后12阶、24阶，以及较显著的36阶、48阶（考虑到模型的简洁性，30阶、40阶可忽略不计）；具有拖尾的特征。前1-11阶不显著因此常规部分可以看成白噪声，而季节部分可以判为AR模型。

（2）模型定阶

从模型识别中我们看到仅季节差分后的序列已经适合建模，但由于增加了一阶差分后序列的趋势分量得到了比较大的改善，而模型复杂度也并未增加太多，因此也可以对一阶差分和季节差分后的序列进行建模。
a) 仅进行季节差分的序列的备选模型
非季节性部分，ACF线性递减，PACF一阶显著，可以定为AR(1)
季节部分，PACF 12阶显著，36阶、48阶较为显著，可以定为SAR(1)或SAR(4)
于是备选模型有:
ARIMA（1，0，0）×（1，1，0）
ARIMA（1，0，0）×（4，1，0）
b) 进行一阶差分和季节差分的序列备选模型
非季节性部分， ACF中滞后1到11阶均不显著，因此可判断为MA(0)。PACF中各滞后阶数也不显著，因此我们将考虑AR(0)；但由于ACF、 PACF在滞后2阶均处于临界状态，在滞后二阶之后相关性下降很快，因此AR的另一个可能的阶数为AR(2)
周期性方面，ACF在滞后12阶之后不再显著，可以假设模型SMA(1)。PACF在滞后4之后不再显著，并具有拖尾特征，所以可以假设模型SAR(4)；由于滞后2阶之后观察到的相关性不像之前的那么大，这意味着它们可能只是来自于噪声，所以也可以考虑模型SAR(2);
所以备选模型有：
ARIMA（0，1，0）×（0，1，1）
ARIMA（0，1，0）×（2，1，0）
ARIMA（0，1，0）×（2，1，1）
ARIMA（0，1，0）×（4，1，0）
ARIMA（0，1，0）×（4，1，1）
ARIMA（2，1，0）×（0，1，1）
ARIMA（2，1，0）×（2，1，0）
ARIMA（2，1，0）×（2，1，1）
ARIMA（2，1，0）×（4，1，0）
ARIMA（2，1，0）×（4，1，1）

（3）参数估计

通过对以上各模型的参数进行估计，得到下面的结果：

模型	参数情况	选择
ARIMA（1，0，0）×（1，1，0）	显著	log likelihood = -443.31, aic = 890.61
ARIMA（1，0，0）×（4，1，0）	显著	log likelihood = -426.36, aic = 862.71
ARIMA（0，1，0）×（0，1，1）	显著	log likelihood = -416.91, aic = 835.83
ARIMA（0，1，0）×（2，1，0）	显著	log likelihood = -437.46, aic = 878.92
ARIMA（0，1，0）×（2，1，1）	不显著
ARIMA（0，1，0）×（4，1，0）	显著	log likelihood = -429.77, aic = 867.54
ARIMA（0，1，0）×（4，1，1）	不显著
ARIMA（2，1，0）×（0，1，1）	显著	log likelihood = -415.02, aic = 836.04
ARIMA（2，1，0）×（2，1，0）	显著	log likelihood = -435.26, aic = 878.52
ARIMA（2，1，0）×（2，1，1）	不显著
ARIMA（2，1，0）×（4，1，0）	显著	log likelihood = -427.5, aic = 867
ARIMA（2，1，0）×（4，1，1）	不显著

根据参数的显著性情况，结合模型的拟合效果和简洁性（对数似然值、AIC值），最终选出2个效果较好的模型，其模型参数如下：

Model1: ARIMA（0，1，0）×（0，1，1）
Coefficients:sma1           -1.0000       s.e.     0.0893   sigma^2 estimated as 1.382:  log likelihood = -416.91,  aic = 835.83  Model2: ARIMA（2，1，0）×（0，1，1）
Coefficients:ar1      ar2     sma1-0.0932   -0.0870   -1.0000s.e.    0.0627   0.0626    0.0924sigma^2 estimated as 1.362:  log likelihood = -415.02,  aic = 836.04

（4）模型诊断

a) 残差分析
Model 1 :

方差齐性：从残差分布图可以看出，残差没有明显趋势，并且几乎所有残差都集中在±3之间，但有一个时间点上的残差超过3接近于4；
正态性：从QQ图可以看到，残差基本符合正态性假设，但尾部有一点较为异常；进行W检验，得到W = 0.99488, p-value = 0.5623，因此可以认为符合正态性假设；
独立性：从ACF图可以看到，在滞后20和28项存在较为显著的相关性；进行游程检验p=0.398，经过Ljung-Box检验p=0.2289，因此在整体上可以认为符合独立性假设。
Model 2：
残差分析图与模型1的结果十分类似不再呈现。于是，与模型1一样，均存在一个异常值点，均能通过W检验、游程检验、Ljung-Box检验，因此符合方差齐性、正态性、独立性假设。
b) 处理异常值
对于残差分析中发现的一个异常值，我们通过软件命令找到，发现它是一个新息异常值，位于第159位（即1986年3月）；于是我们通过对异常值点单独进行回归拟合再加入模型中来消除影响。得到模型3、模型4；模型估计如下：

模型3: 加入IO的ARIMA（0，1，0）×（0，1，1）Coefficients:sma1      IO-1.0000   2.8264s.e.   0.0797   0.8321sigma^2 estimated as 1.322:log likelihood = -411.27, aic = 826.55模型4：加入IO的ARIMA（2，1，0）×（0，1，1）Coefficients:ar1        ar2       sma1      IO-0.0858   -0.0964   -1.0000   2.9400s.e.   0.0628    0.0626     0.0837   0.8623sigma^2 estimated as 1.301:log likelihood = -409.32,aic = 826.65

可以看出与模型1、2相比，对数似然均有所增加，AIC也变小了。
依然进行同样的残差诊断，可以发现：方差齐性、独立性、正态性检验均通过；
进一步检验异常值，此时既没有可加异常值也没有新息异常值。

（5）过拟合

对于模型3的所有尝试得到的过拟合模型，一部分模型并没有有效降低AIC，还有一部分模型降低了AIC但估计的参数值不显著。
对模型4进行过拟合最终获得了表现更好的几个模型；模型ARIMA（2，1，1）×（0，1，1）、模型ARIMA（2，1，2）×（0，1，1）等，比较AIC可以得到ARIMA（2，1，2）×（0，1，1）是最优的，AIC为821.33，对数似然值为-404.67。所以确定这一模型作为最终拟合模型。
模型参数估计如下，记为模型5：
Coefficients:
ar1 ar2 ma1 ma2 sma1 IO
1.7844 -0.9279 -1.8773 1.0000 -1.0000 2.6435
s.e. 0.0250 0.0262 0.0274 0.0289 0.0641 0.8405
sigma^2 estimated as 1.238:log likelihood = -404.67, aic = 821.33
进行最后的残差检验：

根据上面3个图，可以看出：

残差没有明显趋势，并且所有残差都集中在±3之间，符合方差齐性。
残差均分布在QQline附近，基本符合正态性假设；进行W检验，得到W = 0.99534, p-value = 0.6463，因此可以认为符合正态性假设。
ACF图中所有滞后项的相关性均不显著，可以认为残差是白噪声；进行游程检验p= 0.136，经过Ljung-Box检验p= 0.4514，因此可以认为符合独立性假设。
于是，模型5 调整异常值的ARIMA（2，1，2）×（0，1，1）模型通过了所有残差检验。

结果预测

利用第四部分选出的模型5即调整异常值的ARIMA（2，1，2）×（0，1，1）模型，进行预测。由于建模是利用Box-Cox变换过后的数据进行的，所以在画图前需要把模型预测的结果转换成对原序列的预测结果。将原序列的真实值和预测值画在同一图中：

将尾部放大来看，红线表示预测值，黑色虚线表示95%的预测区间，圆圈表示真实值，可以看到，所有的真实值均落在预测值的95%置信区间内，可见这一模型有较高的可信度。

总结

为预测美国新单户住宅的销售量，基于1973年到1995年的月销售量数据建立了乘法季节模型。由于数据集平稳较差，首先采用了Box-Cox变换；转换后的数据显示了一定的趋势和季节性成分，因此进行了季节性的差分、选择性进行了一阶差分。在这之后，通过检查ACF和PACF，根据AIC指标确定了两个候选模型。对于每个模型，进行残差诊断，并去除异常值。继而经过过拟合得到了效果最好的模型：
ARIMA(2,1,2)X(0,1,1)
Coefficients:
ar1 ar2 ma1 ma2 sma1 IO
1.7844 -0.9279 -1.8773 1.0000 -1.0000 2.6435
s.e. 0.0250 0.0262 0.0274 0.0289 0.0641 0.8405
sigma^2 estimated as 1.238:log likelihood = -404.67, aic = 821.33
这一模型的系数显著不为0，有最低的AIC，并且通过了所有的残差诊断，残差为白噪声序列。同时我们看到，它成功地预测了序列的未来数值，真实值在预测值的95%置信区间内。
因此，通过这一模型我们便可以基于过去数据，较为准确地预测未来美国新单户住宅的销量，为各方决策提供参考；进一步根据供求关系，可以帮助我们来判断最佳的购买时期。