时间序列数据可依据变量自身的变化规律,利用外推机制描述序列变化,建立模型并对未来趋势进行预测[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显示,原始序列均值稳定在20左右, 没有明显的上升、下降趋势,序列有较强的季节特征,表现为夏季月份偏高,冬季月份偏低,自相关函数在滞后12阶, 24阶, 36阶…上有强相关性。考虑进行季节差分,消除季节趋势,使用序列变为平稳序列。
>o3_D1=diff(o3, lag=12)#进行季节差分,形成一次季节差分序列
>library(‘forecast’)#载入“forecast”包
>tsdisplay(o3_D1, lag=48, main=‘一次季节差分序列’)#绘制o3_D1序列图、自相关函数图、偏自相关函数图(图 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显示残差的绝大部分自相关函数、偏自相关函数已经落入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年臭氧浓度的预测值和观察值比较见表 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.) |