时间序列ARIMA模型在R中的实现

1、。江M*fJIANGSUUNIVERSITY江苏省餐饮业零售总额分析预测学校：江苏大学学院：财经学院班级：统计1201组员：韩亚琼3120812015马海燕3120812022顾君颖3120812020王培培3120812009陆金龙3120812029白卓3120812028完成时间：2014年12月13日星期六摘要引言数据分析原始数据的获取：本文所有的样本数据均来自江苏统计年鉴一一2014»(得到的样本数据参见表1:表1按行业分社会消费品零售总额年份社会消费品零售总额批发和零售业餐饮业其他行业197884.7979.183.242.37197999.1691.613.903.65

2、1980122.56114.354.723.491981134.79125.165.174.461982150.01138.875.495.651983169.12156.286.146.701984205.05188.807.618.641985262.57240.6910.4511.431986304.58279.2512.5312.801987360.74329.3115.8415.591988471.83432.0620.3319.441989509.56467.1122.4719.981990515.43472.7224.1718.541991578.12529.9427.8620.

3、321992704.52644.6133.6426.271993967.77888.2444.7434.7919941359.611238.3071.4449.8719951741.921573.0195.2173.7019962080.441901.47135.6443.3319972300.612082.71167.9249.9919982453.842208.24192.5253.0819992649.562367.59227.5854.3920002908.462583.19269.5955.6920013233.352845.89326.7160.7620023656.573179.

4、23410.8366.5220034194.503613.67510.9469.8820044892.184333.18496.1021.2220055735.505051.70583.0950.9120066706.195898.79678.8360.8920077985.907023.48810.5669.3820089905.108890.30826.1088.90200911484.1010312.81957.23106.40201013606.8012207.181147.99124.50201115988.3814320.871359.27146.30201218331.30164

5、48.831588.08115.97201320796.5018694.851788.44139.98这里我们仅用到第三列数据，为了方便分析，我们将餐饮业零售总额序列命名为caterts第一步序列的平稳性检验为判断一个序列是否平稳，我们主要通过时序图以及自相关图进行检验。对caterts做时序图，有图形发现有明显的指数趋势，序列非平稳，也可以初步发现江苏省的餐饮业零售总额逐年递增，尤其是在新世纪以后，人们的生活水平逐年提高，对餐饮业的贡献也增大：CATERING1980199020002010gooOTime图1caterts序列时序图因为原序列有明显的指数趋势，故先对数列进行对数变换得到新的

6、数列logcatets,序列图如下，具有明显的非线性增长趋势:图2对数化后的时序图对具有明显线性趋势的数列常用的平稳化措施是差分，我们对logcaterts序列进行一阶差分得到新的数列difflogcaterts,时序图如下：loggedanddiffed图3对数化和一阶差分后的时序图序列进行中心化处理，得到新的数列通过对时序图分析发现数列具有平稳性，为了方便分析，我们对difflogcatertsx。对x进行ADF检验（单位根检验）。R语言中有专门的fUnitRoots包，里面有urdftest功能，是专门进行序列的ADF单位根检验，通过R软件得到如下结果：Title:AugmentedDi

7、ckey-FullerUnitRootTestTestResults:TestregressionnoneCall:lm(formula=乙diffz.lag.1-1+z.diff.lag)Residuals:Min1QMedian3QMax-0.224894-0.0510730.0062610.0432570.242110Coefficients:EstimateStd.ErrortvaluePr(>|t|)z.lag.1-0.592590.19455-3.0460.0047*乙diff.lag0.029090.180670.1610.8731Signif.codes:0*'0

8、.001*'0.01'*'0.05'.'0.1''1Residualstandarderror:0.08821on31degreesoffreedomMultipleR-squared:0.2853,AdjustedR-squared:0.2392F-statistic:6.188on2and31DF,p-value:0.005481Valueoftest-statisticis:-3.046Criticalvaluesforteststatistics:1pct5pct10pcttau1-2.62-1.95-1.61从结果可以看出ADF

9、统计量为-3.046,在1%,5%,10%的置信水平下均拒绝原假设，认为序列x平稳。第二步模型的识别与定阶模型的识别与定阶主要是通过对序列的自相关和偏自相关观察得到的。通过R软件的到序列x的自相关和偏自相关图：图4序列x的自相关图Seriesx图5序列x的偏自相关图通过自相关图和偏自相关图，可以发现序列x自相关明显拖尾，偏自相关有结尾的性质，在第十阶的时候，PACF超出了两倍的标准差，这可能是由于偶然因素引起的，可以考虑使用AR(1)和AR(11)模型，但是为了数据的完整性，尽量避免滞后期，可以考虑AR(1)模型。在R语言中有专门的forecast包，包里的auto.arima函数可以对序列模

10、型进行智能识别建模。通过R语言的智能识别功能得到如下结果：Series:xARIMA(1,0,0)withzeromeanCoefficients:ar10.4126s.e.0.1512sigmaA2estimatedas0.006901:loglikelihood=37.33AIC=-70.65AICc=-70.28BIC=-67.54通过如剑的输出结果我们发现结果与我们预测的一样，均是AR(1)模型。得到的模型方程为：?=0.1426?一1+?第三步模型的适应性检验模型的适应性检验实质上就是检验残差序列?是否为白噪声序列。其中最主要的是?外序列的独立性检验。通过R软件得到模型的残差序列的时

11、序图(偏)自相关图：Residuals图6残差序列的时序、自相关和偏自相关图通过对图形的观察，初步判断残差序列?分为白噪声。R语言中的Box.test函数也可对残差序列进行白噪声检验。且得到的结果如下:LAGLBp1,10.021283920.88400812,20.084847640.95846353,30.437509620.93238944,40.696146050.95180335,52.029133050.84509946,62.032049590.91672757,72.248820310.94479578,82.454730560.96380859,98.500170790.48

12、4629310,1010.925501510.363351511,1110.939876320.448313912,1213.844479050.3107497可以发现无论滞后几期，P值均大于0.05,所以接受原假设，认为残差序列通过纯随机性检验。故模型最终为ARIMA(1,0,0),即为AR(1)第四步模型参数的估计模型参数的估计主要有矩估计，最小二乘估计（ML估计），极大似然估,这里我们主要用的是最小二乘估计，估计的结果在第二步模型的识别中R语言智能建模中给出了，这里不做重复论述。模型的估计方程如下：?=0.1426?一1+?第五步模型的预测将模型的拟合值与模型的真实值进行画图对比，得到如

13、下图:19SC198519901995200020052010Time图7拟合值一真实值发现拟合值与真实值之间有一定差别,这是由于残差序列的存在，将拟合值与残差序列相加得到的结果与X序列作图，结果如下：图8加上残差的拟合值一真实值观察图8,可以发现拟合的特别好，所有的点都落在了真实值上。我们用拟合的有效模型进行短期预测，比如我们预测2014年、2015年、2016、2017年和2018年的餐饮业零售总额。先预测2014年、2015年、2016年、2017年和2018年的x,再预测餐饮业零售总额。在R语言中可以使用predict函数对序列进行预测，这里我们进行向前5步预测,得到点预测，区间预测。

14、结果如下：$predTimeSeries:Start=2014End=2018Frequency=11-0.0254059516-0.0104836538-0.0043260335-0.0017851186-0.0007366213$seTimeSeries:Start=2014End=2018Frequency=110.083071630.089866320.090972740.091159800.09119162利用公式?土？告XS.E.2可以得到预测值的区间估计，点估计和区间估计具体结果如下:序列X20142015201620172018点预测值-0.0254060-0.0104837-

15、0.004326-0.0017851-0.0007366下界值(L)-0.1882263-0.1866216-0.1826326-0.1804583-0.1794722上界值(U)0.13741440.16565430.17398050.17688810.1779990通过还原零均值、差分和对数变换，就可以得到江苏省餐饮业零售总额(caterts(2014至U2015年的预测值：利用公式：?1?+?jrIoLjnIJ1IJ(bbk?1?exp7?0?得到江苏省2014年餐饮业的零售手总额预测值，以此类推，分别得到2015,2016,2017,2018年餐饮业零售总额的预测值。值得一体的是时间序

16、列预测效果越好，则向前预测的阶数就越少。所以在该模型中预测最好的是2014年。caterts20142015201620172018点预测值2088.2452474.9672951.4233528.5564222.97值得一体的是时间序列预测效果越好，则向前预测的阶数就越少。所以在该模型中预测最好的是2014年。4、 参考文献1王振龙.应用时间序列分析M.北京：中国统计出版社，20102 PaulTeetor.R语言经典实例M.北京：机械工业出版社，20133 RobertI.Kabacoff.R语言实战M.北京：人民邮电出版社，20135、 附录R语言代码：'''rc

17、atering<-read.table("E:catering.txt",header=F)caterts<-ts(catering,start=1978,frequency=1)#数据时间序列化plot.ts(caterts,main=("CATERING"),ylab="CATERING")#序列图logcaterts<-log(caterts)difflogcaterts<-diff(log(caterts)#对数差分plot(difflogcaterts,main="loggedanddiff

18、ed")x<-scale(difflogcaterts,T,F)#零均值化acf(x,10)#自相关pacf(x,10)#偏自相关corr<-acf(x,10)cov<-acf(x,10,type="covariance")library(fUnitRoots)#ADF单位根检验urdfTest(x)library(forecast)#建模fit<-auto.arima(x)fitresid<-fit$residual#生成残差序列residplotForecastErrors(fit$residuals)# 白噪声检验Result=

19、0LAG=0LB=0p=0for(iin1:12)Btest=Box.test(resid,type="Ljung-Box",lag=i)LAGi=iLBi=Btest$statisticpi=Btest$p.valueResult=cbind(LAG,LB,p)Result#接受原假设，认为resid序列为白噪声# 预测pre=predict(fit,n.ahead=5)preU=pre$pred+1.96*pre$seL=pre$pred-1.96*pre$sets.plot(x,pre$pred,col=1:2)lines(U,col="blue"

20、,lty="dashed")lines(L,col="blue",lty="dashed")plot(x,type="o")lines(fitted(fit),col="2",pch=7)points(fitted(fit),col="2",pch=7)plot(x,type="o")points(fitted(fit)+resid,col="blue",pch=20)# 预测20142015年x2014<-0.0254059516+logcaterts36+mean(difflogcaterts)x2015<-0.0104836538+x2014+mean(difflogcaterts)x2016<-0.0043260335+x2015+mean(difflogcaterts)x2017<-0.0017851186+x2016+mean(difflogcaterts)x2018<-0.0007366213+x2017+mean(difflogcaterts)forecast1<-exp(x2014