时间序列数据可依据变量自身的变化规律，利用外推机制描述序列变化，建立模型并对未来趋势进行预测[1]。ARIMA乘法季节模型是研究这种时间序列数据的经典方法[2]，R软件作为一种自由、开源和免费的软件，应用越来越广泛，利用R软件实现ARIMA乘法季节模型的报道多为应用性研究, 缺乏详细操作步骤和流程，本文以大气污染物臭氧(O3)浓度数据为例, 就ARIMA乘法季节模型的R软件实现做一介绍。
1 ARIMA乘法季节模型简介ARIMA乘法季节模型是ARIMA模型的一种特殊形式。ARIMA模型全称为自回归移动平均混合模型(Autoregressive Integrated Moving Average Models)。包括自回归模型(autoregressive model, AR)、差分(Integrated, I)、移动平均模型(moving average model, MA)3个部分，是美国统计学家Geogre E. P. Box和英国统计学家GunlymM. Jenkins创立的时间序列分析方法中的基本模型之一[2]。
1.1 AR、MA、ARMA(p, q)、ARIMA(p, d, q)、ARIMA(p, d, q)×(P, D, Q)s模型[2-3]对于时间序列数据xt，AR模型具有如下结构
| ${{\rm{x}}_{\rm{t}}} = {{\rm{ \mathsf{ φ} }}_{\rm{1}}}{{\rm{x}}_{{\rm{t - 1}}}} + {{\rm{ \mathsf{ φ} }}_{\rm{2}}}{{\rm{x}}_{{\rm{t - 2}}}} + \cdots + {{\rm{ \mathsf{ φ} }}_{{\rm{pxt - p}}}} + {{\rm{u}}_{\rm{t}}} $ |
该模型表示，t时刻的观察值与前p个时刻的观察值呈线性关系, 称为p阶自回归，φ为自回归系数，记为AR(p)。ut表示残差序列，又称白噪声(white noise)序列。用后移动算子B表示为φ(B)xt=ut，其中φ(B)=1-φ1B-φ2B2-…-φpBp；MA模型具有如下结构
| ${{\rm{x}}_{\rm{t}}} = {{\rm{u}}_{\rm{t}}} - {{\rm{ \mathsf{ θ} }}_{\rm{1}}}{{\rm{u}}_{{\rm{t - 1}}}} - {{\rm{ \mathsf{ θ} }}_{\rm{2}}}{{\rm{u}}_{{\rm{t - 2}}}} - \cdots - {{\rm{ \mathsf{ θ} }}_{{\rm{qut - q}}}}, $ |
表示xt在t时刻的取值与q个随机扰动呈线性关系，称为q阶移动平均模型，θ为滑动平均系数，记为MA(q)，用后移动算子B表示为xt=θ(B)ut，其中θ(B)=1-θ1B-θ2B2-…-θqBq；对于AR(p)模型，若ut不是白噪声，表现为MA(q)的形式, 则为ARMA(p, q)模型，即：φ(B)xt=θ(B)ut。ARMA(p, q)模型要求做预测的时间序列为平稳序列，图形显示所有观察值围绕某一水平直线上下随机波动，若序列不平稳，需经过差分变换，使其平稳，即ARIMA(p, d, q)模型，d为差分次数。
某些以月和季度为单位的数据，如大气污染物的月均浓度，常具有明显的季节特征，季节模型的时间单位为相应的周期S, 表示为ARIMA(P, D, Q)s，即Φ(BS)Wt=Θ(BS)ut，P、D、Q分别表示以s为间距的自回归、差分、移动平均阶数，令Wt=▽d▽SDxt，表示经过差分和季节差分得到的序列, 其中Φ(BS)=1-Φ1BS-Φ2B2S-…-ΦPBPS，Θ(BS)=1-Θ1BS-Θ2B2S-…-ΘQBQS，上述全部模型综合在一起即为ARIMA乘法季节模型，表示为：
| ${{\rm{ \mathsf{ φ} }}_{\rm{p}}}(\rm{B}){{\rm{\Phi }}_{\rm{p}}}({{\rm{B}}^{\rm{S}}}){\rm{W}_t} = {{\rm{ \mathsf{ θ} }}_{\rm{q}}}({\rm{B}}){{\rm{\Theta }}_{\rm{Q}}}({B^S}){u_t} $ |
记为: ARIMA(p, d, q)×(P, D, Q)s，模型需要对p、d、q、P、D、Q进行定阶，从而进行估计。
1.2 模型的分析步骤[2, 4] 1.2.1 序列平稳化和纯随机性检验通过序列图、自相关系数函数图(auto correlation function，ACF)、偏自相关系数函数图(partial autocorrelation function, PACF)等分析序列平稳性, 查看是否有趋势性、季节性等特征，并对序列进行ADF检验，考察平稳性。若序列不平稳，需要进行差分和(或)季节差分，将序列转换成平稳序列；利用LB(Ljung-Box)统计量对序列进行白噪声检验：若为非白噪声，说明序列间存在相关关系，可以进行分析。
1.2.2 模型识别与定阶根据差分次数、自相关函数、偏自相关函数等初步确定p、d、q、P、D、Q的阶数，基本原则是：ACF图拖尾，PACF图p阶截尾，选择AR(p)模型；ACF图q截尾，PACF图拖尾，选择MA(q)模型；ACF图、PACF图均拖尾选择ARMA(p, q)模型。
1.2.3 参数估计与模型诊断利用极大似然估计、最小二乘估计等方法，估计模型参数，参数估计后，对模型的残差进行白噪声检验，判断模型的适用性。依据赤池信息准则(AIC)、贝叶斯信息准则(BIC)等方法调整模型阶数，对模型进行比较和优化，选取较为简洁的最佳模型。
1.2.4 预测应用用选定的模型对序列进行预测。
2 材料与方法利用美国芝加哥市大气污染物臭氧(O3)数据为例，进行模拟，数据来源于R自带的程序包，数据库汇总了1987—2000年芝加哥市逐日O3浓度，将其整理成月均浓度数据分析(表 1)。其中，1987年1月—1999年12月的数据用于建立模型，所建模型用来预测2000年(1—12)月的月均浓度，将预测值和实际观察值进行比较，计算相对误差，对模型进行评估。
| 年份/年 | 月份/月 | 臭氧浓度/(μg/m3) |
| 1987 | 1 | 7.7 |
| 1987 | 2 | 11.9 |
| 1987 | 3 | 23.4 |
| … | … | … |
| 2000 | 10 | 16.4 |
| 2000 | 11 | 8 |
| 2000 | 12 | 10.4 |
3 R软件实现[2, 5] 3.1 序列平稳性和纯随机性检验
>chi < -read.csv("chi.csv", header=TRUE, sep=', ') #导入数据
>o3 < -ts(chi＄o3, start=c(1987, 1), frequency=12)#将臭氧浓度转为时间序列变量，频率为12，即为月度数据。
>plot(o3, type=‘o’, ylab=“臭氧浓度”, xlab=“年份”, main=“原始序列图”) #绘制序列图
>acf(as.vector(o3), lag=36, main=“原始序列”)#绘制序列自相关图(图 1)。
|
| 图 1 原始序列图 |
图 1显示，原始序列均值稳定在20左右, 没有明显的上升、下降趋势，序列有较强的季节特征，表现为夏季月份偏高，冬季月份偏低，自相关函数在滞后12阶, 24阶, 36阶…上有强相关性。考虑进行季节差分，消除季节趋势，使用序列变为平稳序列。
>o3_D1=diff(o3, lag=12)#进行季节差分，形成一次季节差分序列
>library(‘forecast’)#载入“forecast”包
>tsdisplay(o3_D1, lag=48, main=‘一次季节差分序列’)#绘制o3_D1序列图、自相关函数图、偏自相关函数图(图 2)。
|
| 图 2 一次季节差分序列 |
图 2序列图、自相关函数图显示经过季节差分，原始序列明显的季节趋势已经基本消除，可以认为o3_D1序列为平稳序列。
>library(tseries)#载入“tseries”程序包
>adf.test(o3_D1)#单位根检验
P＜0.01，拒绝序列为非平稳序列的原假设，可以认为一次季节差分后的序列为平稳序列。
>Box.test(o3, lag=6, type="Ljung-Box")；Box.test(o3, lag=12, type="Ljung-Box")#纯随机性检验
LB统计量结果结果显示，滞后6阶和滞后12阶，P值均小于0.01，在α=0.05水平，拒绝序列为纯随机的原假设，可以认为原始序列为非白噪声序列，可以进行后续分析
3.2 模型定阶根据图 2，季节差分后的序列，ACF图在1阶、12阶后，PACF图在2阶、24阶后快速衰减，表现出截尾特征，初步选定p=2、d=0、q=1、P=2、D=1、Q=1。
3.3 模型拟合与优化>fit1 < -arima(o3, order=c(2, 0, 1), seasonal=list(order=c(2, 1, 1), period=12), method=‘ML’) ##拟合模型, 选用极大似然估计的方法
>acf(as.vector(fit1＄residuals), lag=36, main=“残差序列”);
pacf(as.vector(fit1＄residuals), lag=36)#绘制残差的ACF、PACF图(图 3)。
|
| 图 3 残差序列 |
图 3显示残差的绝大部分自相关函数、偏自相关函数已经落入2倍标准差之内，可以认为残差序列为白噪声。
>Box.test(fit1＄residuals, lag=6, type="Ljung-Box"); Box.test(fit1＄residuals, lag=12, type="Ljung-Box")#残差的白噪声检验
LB统计量结果结果显示，滞后6阶和滞后12阶，P值均大于0.05，在α=0.05水平，可以认为残差序列为白噪声序列, 模型1显著有效。
>fit2 < -arima(o3, order=c(0, 0, 1), seasonal=list(order=c(0, 1, 1), period=12), method=‘ML’)#经过多次试探性降阶尝试，拟合AIC值较小、更简洁的模型2
AIC(fit1, fit2)#比较两个模型的AIC值
结果显示AIC2(701.9)＜AIC1(707.3)，选择模型2作为最终模型。
3.4 模型预测>o3.fore=forecast(fit2, h=12, level=0.95)#利用模型2进行一个周期，即2000年预测
>summary(o3.fore)#查看预测结果
>plot(o3.fore)#绘制预测图(图 4)。
|
| 图 4 2000年预测图 |
4 结果
2000年臭氧浓度的预测值和观察值比较见表 2，结果显示平均相对误差为5.6%。
| 年份 | 月份 | 观测值 | 预测值 | 绝对误差 | 相对误差/% |
| 2000 | 1 | 12.8 | 12.3 | -0.5 | -3.9 |
| 2000 | 2 | 17.6 | 16 | -1.6 | -9.1 |
| 2000 | 3 | 20.9 | 23.2 | 2.3 | 11 |
| 2000 | 4 | 21.4 | 23.1 | 1.7 | 7.9 |
| 2000 | 5 | 23.9 | 26.9 | 3 | 12.6 |
| 2000 | 6 | 25.3 | 30.1 | 4.8 | 19 |
| 2000 | 7 | 28.2 | 28.9 | 0.7 | 2.5 |
| 2000 | 8 | 25.3 | 23.8 | -1.5 | -5.9 |
| 2000 | 9 | 19.4 | 20.8 | 1.4 | 7.2 |
| 2000 | 10 | 16.4 | 13.2 | -3.2 | -19.5 |
| 2000 | 11 | 8 | 11.9 | 3.9 | 48.8 |
| 2000 | 12 | 10.4 | 10 | -0.4 | -3.8 |
| 平均相对误差 | 5.6 |
5 小结
本文对ARIMA乘法季节模型做了简要的介绍，以有季节特征的臭氧浓度数据为例，给出了用R软件实现模型拟合、检验、图形绘制以及预测的详细过程。R软件由于其较大的灵活性，近年来在数据分析、作图等领域应用非常广泛，不同的统计方法通常都能找到相应的软件包，便于使用者使用，本文使用了forecast、tseries两个软件包。ARIMA模型通过调整相关参数，也可以实现以日、周为单位的数据分析[6]。除大气污染物的数据外，类似具有相关特征的数据也可采用此方法进行预测，如医院的门诊量[7]、日死亡数、传染病的发病数据[8]等，在实际工作中使用者应根据数据特征进行合理分析。
| [1] |
温亮, 张秀山, 李承毅, 等. 季节分解法和ARIMA法预测乌鲁木齐市肺结核发病趋势效果分析[J]. 军事医学, 2017, 41(4): 287-290. (In English: Wen L, Zhang XS, Li CY, et al. Seasonal decomposition and ARIMA methods in prediction of tuberculosis incidence in Urumqi, China[J]. Mil Med Sci, 2017, 41(4): 287-290.) |
| [2] |
王燕. 应用时间序列分析[M]. 3版. 北京: 中国人民大学出版社, 2012. (In English: Wang Y. Applied Time Series Analysis[M]. 3rd ed. Beijing: China Renmin University Press, 2012.)
|
| [3] |
张文彤, 董伟. SPSS统计分析高级教程[M]. 2版. 北京: 高等教育出版社, 2013.
|
| [4] |
胡跃华, 廖家强, 冯国双, 等. 自回归移动平均模型在全国手足口病疫情预测中的应用[J]. 疾病监测, 2014, 29(10): 827-832. (In English: Hu YH, Liao JQ, Feng GS, et al. Application of multiple seasonal autoregressive integrated moving average model in prediction of incidence of hand foot and mouth disease in China[J]. Dis Surve, 2014, 29(10): 827-832. DOI:10.3784/j.issn.1003-9961.2014.10.018) |
| [5] |
Cryer JD, Chan KS. Time Series Analysis with Applications in R[M]. Beijing: China Machine Press, 2011.
|
| [6] |
牟敬锋, 赵星, 樊静洁, 等. 基于ARIMA模型的深圳市空气质量指数时间序列预测研究[J]. 环境卫生学杂志, 2017, 7(2): 102-107, 117. (In English: Mou JF, Zhao X, Fan JJ, et al. Time series prediction of AQI in Shenzhen based on ARIMA model[J]. J Environ Hyg, 2017, 7(2): 102-107, 117.) |
| [7] |
杨帆, 秦银河, 刘丽华. ARIMA模型在门诊人次预测中的应用[J]. 中华医院管理杂志, 2009, 25(1): 28-31. (In English: Yang F, Qin YH, Liu LH. Application of ARMA Model in prediction of outpatient headcount[J]. Chin J Hosp Admin, 2009, 25(1): 28-31. DOI:10.3760/cma.j.issn.1000-6672.2009.01.011) |
| [8] |
王建书, 刘强, 覃江纯, 等. 基于ARIMA乘积季节模型的苏州市介水传染病发病预测研究[J]. 环境卫生学杂志, 2017, 7(6): 417-420. (In English: Wang JS, Liu Q, Qin JC, et al. Prediction of incidence for water-borne diseases on a multiple seasonal ARIMA model in Suzhou[J]. J Environ Hyg, 2017, 7(6): 417-420.) |