1、第10章 时间序列分析2020/4/10常用时间序列模型ARIMA 模型时间序列对象时间序列对象在R软件中，使用时间序列建模前需要先将数据存储到一个时间序列对象中。我们可以使用函数ts()将数值类型的观测对象存储为时间序列对象。使用格式：ts(data=NA,start=1,end=numeric(),frequency=1,)其中，data是时间序列观测值对象，必须为数值类型的向量、矩阵或数据框；start是用来指定时间序列观测值对象的第一个时间点，比如2000年1月，则设置start=c(2000,1)；end用来指定时间序列的终止时间点；frequency用来指定数据在一年中的频数。还可
2、以通过函数as.ts()可以将对象转换成时间序列；通过函数is.ts()可以判断对象是否为时间序列对象。时间销量/元0510152025303530004000ARIMA 模型绘制时间序列图绘制时间序列图R软件中，可以使用plot.ts()函数来画出时间序列的时序图。plot.ts()用法同plot。根据平稳时间序列的均值和方差都为常数的性质，平稳序列的时序图显示该序列值始终在一个常数附近随机波动，而且波动的范围有界；如果有明显的趋势性或者周期性那它通常不是平稳序列。ARIMA 模型时间序列检验分析（1）自相关性检验自相关图中的两条虚线标示置信界限是自相关系数的上下界。如果自相关系迅速衰减落入
3、置信区间内，就可能是白噪声；如果自相关系数超出置信区间，那么表示存在相关关系，而且从哪一阶落在置信区间内，就表示自相关的阶数是几阶。使用格式：acf(x,lag.max=NULL,type=c(correlation,covariance,partial),plot=TRUE,na.action=na.fail,demean=TRUE,)pacf(x,lag.max,plot,na.action,)acf()函数中，参数 x为观测值序列，acf()为观测值序列自相关函数，lag.max为与acf对应的最大延迟，type为计算acf的形式，默认为 correlation。当没有输出，即为 acf
4、(Series)时，画观测值序列的自相关系数图。pacf()函数中的输入参数与输出参数的含义同 acf()函数的类似。在 acf()和pacf()中设定plot=FALSE可以得到自相关和偏自相关的真实值。ARIMA模型原始序列时序图原始序列的自相关图时序图时序图显示该序列具有明显的单调递增趋势，可以判断为是非平稳序列；自相关图自相关图 显示自相关系数长期大于零，说明序列间具有很强的长期相关性，可以判断为非平稳序列；ARIMA 模型时间序列检验分析（2）单位根检验对时间序列的平稳性检验通常使用单位根检验的方法。在R软件中，单位根检验使用fUnitRoots包中的unitrootTest()函数
5、可以实现使用格式：unitrootTest(x,lags=1,type=c(c,c,ct),title=NULL,description=NULL)其中，输入参数x为观测值序列，lags为用于校正误差项的最大滞后项，type为单位根的回归类型，返回的参数p值，p值小于0.05表示满足单位根检验。单位根检验统计量对应的 p值显著大于 0.05，判断该序列为非平稳序列（非平稳序列一定不是白噪声序列）。ARIMA 模型ARIMA建模分析（1）非平稳时间序列差分对于非平稳时间序列，首先需要对其进行差分直到得到一个平稳时间序列。在R软件中，可以使用diff()函数对时间序列进行差分运算。diff()函数
6、的使用格式：diff(x,lag=1,differences=1,.)其中，输入参数硜 代表观测值序列；汜条代表差分运算的步数，缺省值代表一步差分；?楤晦牥湥散屳代表差分运算的阶数，缺省值代表一阶差分。对一阶差分后的序列再次做平稳性判断过程同上。时间销量残差/元5101520253035-200-1000100200300051015202530-0.20.00.20.4Series difsalesLagACF一阶差分之后序列的时序图一阶差分之后序列的自相关图结果显示，一阶差分之后的序列的时序图在均值附近比较平稳的波动、自相关图有很强的短期相关性、单位根检验p值小于0.05，所以一阶差分之后
7、的序列是平稳序列。ARIMA模型（2）时间序列模型识别定阶使用 R软件中的acf()和pacf()函数来分别给出时间序列的自相关图和偏自相关图。可根据自相关图和偏自相关图对时间序列模型进行定阶。根据自相关图和偏自相关图对时间序列模型进行定阶：?若平稳序列的偏相关系数是截尾的，而自相关系数是拖尾的，则序列适合AR模型；?若平稳序列的偏相关系数是拖尾的，而自相关系数是截尾的，则序列适合MA模型；?若平稳序列的偏相关系数与自相关系数都是拖尾的，则序列适合ARMA模型；ARIMA 模型051015202530-0.20.00.20.4LagPartial ACFSeries difsales05101
8、5202530-0.20.00.20.4Series difsalesLagACFpacf(difsales,lag.max=30)一阶差分后序列的偏自相关图一阶差分后自相关图中，ACF值在一阶后迅速跌入置信区间，并且数值徘徊在置信区间，没有收敛趋势，显示出拖尾性。偏自相关图中，PACF值在一阶后迅速跌入置信区间，并且有向零收敛的趋势，显示出截尾性，所以可以考虑用AR模型拟合1阶差分后的序列，即对原始序列建立 ARIMA(1,1,0)模型。一阶差分之后序列的自相关图acf(difsales,lag.max=30)ARIMA模型另外，模型还可以通过 BIC进行定阶。此处计算 ARMA(p,q)当
9、p和q均小于等于 5的所有组合的BIC信息量，取其中 BIC信息量达到最小的模型阶数。BIC(Intercept)test-lag1test-lag2test-lag3test-lag4test-lag5error-lag1error-lag2error-lag3error-lag4error-lag5171513129.37.65.73.5一阶差分后序列的 BIC图显示BIC值从下往上，依次递减。模型选用变量的单元格用阴影表示。较好的模型（具有较低的 BIC值）处于较高的行中。第一行中，test-lag1被选入模型，error-lag1到error-lag5均未被选取，取零阶。即在 test
10、-lag1，error-lag0处，取得最小 BIC值。因此p值为1、q值为0时，BIC值最小。p、q定阶完成,，即对原始序列建立 ARIMA(1,1,0)模型。ARIMA 模型（3）ARIMA模型R软件中，arima()函数设置时序模式的建模参数，创建ARIMA时序模型或者把一个回归时序模型转换为ARIMA模型。使用格式：Arima(Series,order,seasonal,period,method,.)其中，Series为观测值序列，order为构建的ARIMA(p,d,q)模型的参数，seasonal为模型的季节性参数，period为观测值序列的周期，method为估计模型参数所使用
11、的方法。ARIMA 模型（4）白噪声检验在使用ARIMA模型定阶完成后，还要对模型进行假设检验，检验残差序列是否为白噪声序列。R软件中，Box.test()函数用于检测序列是否符合白噪声检验。使用格式：Box.test(x,lag=1,type=c(Box-Pierce,Ljung-Box),fitdf=0)其中，x为Arima()函数返回的结果对象的residuals残差。ARIMA 模型p-value 0.05，说明残差为白噪声序列，通过白噪声检验。ARIMA 模型预测（1）模型预测根据参数估计得到的时间序列模型对未来的序列值进行预测，通过forecast包中的函数forecast.Ari
12、ma()来完成。使用格式：forecast.Arima(object,h=10,level=c(80,95),)object 是函数arima()返回的对象；h指定预测的时间点；level指定预测区间的置信水平，默认情况为80%和95%置信水平下的预测区间。Forecast的输入结果是指：通过 arima模型，置信水平为 80%和95%时，预测出未来五天的营业额及营业额区间。ARIMA 模型预测（2）绘制原始值与预测值图形为查看原始值与预测值的图形效果，可以使用plot.forecast()函数绘制原始值与预测值的图形。使用格式：plot.forecast(x,)其中，x为forecast.A
13、rima()函数返回的结果对象。Forecasts from ARIMA(1,1,0)010203040300040005000使用函数plot.forecast()可以查看原始值和预测值的图形效果，例如从图中可以直观地看到营业额的发展趋势。plot.forecast(forecast)其他时间序列模型（1）组合模型时间序列的变化主要受到长期趋势、季节变动、周期变动和噪声变动这四个因素的影响。根据序列的特点，可以构建加法模型和乘法模型。decompose()函数、stl()函数可以估计出时间序列中趋势的、季节性的和不规则的部分，而此时间序列需是可以用相加模型描述的。函数使用格式：decompo
14、se(x,type=c(additive,multiplicative),filter=NULL)在decompose()函数中，x时间序列对象；type指定分解为加法模型还是乘法模型；filter是滤波系数。stl(x,s.window,s.degree=0,)在stl()函数中，x同样为时间序列对；s.window因为没有默认值，所以必须手动设置，可以采用 periodic或Loess方法提取季节跨度，若采用Loess方法，s.window的值必须为大于7的奇数；s.degree可取1或0，为局部多项式拟合季节性提取的程度。sales sales.de=decompose(sales,ty
15、pe=additive)sales.de$xJan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec2013 3023 3039 3056 3138 3188 3224 3226 3029 2859 2870 2910 30122014 3142 3252 3342 3365 3339 3345 3421 3443 3428 3554 3615 36462015 3614 3574 3635 3738 3707 3827 4039 4210 4493 4560 4637 47552016 4817$seasonalJan Feb Mar Apr May J
16、un Jul Aug Sep Oct Nov Dec2013 40 33 50 44-56-65 79 22-94-50-23 202014 40 33 50 44-56-65 79 22-94-50-23 202015 40 33 50 44-56-65 79 22-94-50-23 202016 40$trendJan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec2013 NA NA NA NA NA NA 3053 3067 3087 3109 3125 31362014 3149 3174 3215 3268 3325 3381 3427 34
17、60 3486 3514 3545 35802015 3626 3684 3760 3846 3931 4020 4116 NA NA NA NA NA2016 NA$randomJan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec2013 NA NA NA NA NA NA 94 -59-134-189-192-1442014 -47 44 77 54 69 29 -86 -39 36 90 93 462015 -52-143-175-152-168-127-156 NA NA NA NA NA2016 NA$figure1 40 33 50 44-
18、56-65 79 22-94-50-23 20$type?慜摤瑩癩履attr(,class)?摜捥浯潰敳?獴30004500observed32004000trend-1000seasonal-200-501002013.02013.52014.02014.52015.02015.52016.0randomTimeDecomposition of additive time series plot(sales.de)图中包含了四部分，自上而下依次为：原始时间序列的观测值、时间序列分解趋势图、时间序列分解季节变动图、时间序列分解噪声图。其他时间序列模型30004500data-40040seas
19、onal35004500trend-20002013.02013.52014.02014.52015.02015.52016.0remaindertime sales.stl=stl(sales,s.window=periodic)sales.stlCall:stl(x=sales,s.window=periodic)Componentsseasonal trend remainderJan 2013 5.73 3109 -91.9Feb 2013 -53.40 3103 -10.3Mar 2013 -28.08 3096 -12.2Apr 2013 0.99 3091 45.7May 201
20、3 -41.61 3086 143.4Jun 2013 -26.07 3085 165.6 Mar 2015 -28.08 3752 -89.2Apr 2015 0.99 3843 -105.6May 2015 -41.61 3933 -184.3Jun 2015 -26.07 4036 -182.5Jul 2015 32.13 4138 -131.4Aug 2015 -4.78 4243 -28.0Sep 2015 -7.68 4347 153.5Oct 2015 19.36 4456 85.1Nov 2015 37.73 4564 35.4Dec 2015 65.68 4675 13.8J
21、an 2016 5.73 4787 24.2 plot(sales.stl)stl季节性时间序列分解图图中包含了四部分，自上而下依次为：原始时间序列的观测值、时间序列分解趋势图、时间序列分解季节变动图、残差自相关图。采用decompose()与stl()对时间序列进行分解后，建模、预测的工作一般通过指数平滑算法实现。其他时间序列模型（2）指数平滑法在R软件中，简单指数平滑法、Holt双参数线性指数平滑法和Winters线性和季节性指数平滑法可以通过HoltWinters()函数实现。HoltWinters()函数会给出时间序列自身时间段内的预测，生成的预测值存在一个变量列表fitted中。可以
22、用plot()函数绘制原始数据和预测值的图形查看预测效果。使用格式：HoltWinters(x,alpha=NULL,beta=NULL,gamma=NULL,seasonal=c(additive,multiplicative),start.periods=2,l.start=NULL,b.start=NULL,s.start=NULL,optim.start=c(alpha=0.3,beta=0.1,gamma=0.1),optim.control=list()其中，x为时间序列对象；alpha、beta、gama都是HoltWinters算法的过滤系数。beta和gama是Holt指数平
23、滑法或Winters指数平滑法的参数。如果beta设置为FALSE，该函数将做指数平滑；如果gama指定为FALSE，那么该函数适用于拟合非季节性模型；如果gama与beta同时指定为FALSE，那么该函数适用于简单指数平滑模型；通过seasonal参数选择”摡楤楴敶 或浜汵楴汰捩瑡癩履季节性模型，默认选择”additive”而且仅在gamma非零时生效；start.period是用于x对象的frequency自动检测，不能小于2；l.start、b.start、s.start分别表示启动值、趋势值和季节分量的初始值；optim.start设置向量命名的组件alpha、beta、gamma包含
24、优化的初始值，必须指定唯一需要的值，在只使用alpha、beta、gamma中的一个参数时忽略本参数。HoltWinters()函数仅得到预测模型，如果要对未来的时间做预测，需要调用forecast()函数。其他时间序列模型 hw.sales=HoltWinters(sales,alpha=TRUE,beta=TRUE,gamma=TRUE)hw.salesHolt-Winters exponential smoothing with trend and additive seasonal component.Call:HoltWinters(x=sales,alpha=TRUE,beta=T
25、RUE,gamma=TRUE)Smoothing parameters:alpha:TRUEbeta:TRUEgamma:TRUECoefficients:,1a 4790.861111b -54.833333s1 110.763889s2 159.805556s3 130.597222s4 46.722222s5 -3.069444s6 206.388889s7 -4.444444s8 -195.236111s9 -205.611111s10-181.361111s11 -90.694444s12 26.138889 plot(hw.sales)Holt-Winters filteringT
26、imeObserved/Fitted2014.02014.52015.02015.52016.0300040005000餐饮营业额Winters模型拟合图该餐馆的营业数据存在明显的趋势性与季节性。因此可以采用Winters线性和季节性指数平滑法进行建模预测。其他时间序列模型0.00.20.40.60.81.0-0.40.00.40.8LagACFSeries hw.model$residuals建模并对模型的残差进行自相关检验与白噪声检验。library(forecast)hw.model=recast.HoltWinters(hw.sales,h=6,level=c(80,95)hw.mod
27、el Point Forecast Lo 80 Hi 80 Lo 95 Hi 95Feb 2016 4846.792 4640.947 5052.636 4531.980 5161.603Mar 2016 4841.000 4380.718 5301.282 4137.060 5544.940Apr 2016 4756.958 3986.759 5527.157 3579.041 5934.876May 2016 4618.250 3490.794 5745.706 2893.955 6342.545Jun 2016 4513.625 2987.043 6040.207 2178.919 68
28、48.331Jul 2016 4668.250 2704.620 6631.880 1665.138 7671.362acf(hw.model$residuals)Box.test(hw.model$residuals,lag=10,type=Ljung-Box)Box-Ljung testdata:hw.model$residualsX-squared=4.3666,df=10,p-value=0.9293Winters模型残差自相关图由Winters模型残差自相关图可知，各阶的残差系数都没有超过执行区间，可以初步判断该残差不是自相关的；通过白噪声检验P值为0.9293，不能拒绝原假设，说明残差为白噪声序列，不存在自相关性。hw.model模型对该餐饮营业额拟合的效果优秀。其他时间序列模型通过模型残差的检验后，画出原始值与预测值的图形，以便直观观察未来半年的餐饮营业额的变化趋势。plot(hw.model)Forecasts from HoltWinters2013.02014.02015.02016.0200040006000餐饮营业额 Winters模型预测图其他时间序列模型
