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

《时间序列ARIMA模型在R中的实现》由会员分享，可在线阅读，更多相关《时间序列ARIMA模型在R中的实现（11页珍藏版）》请在装配图网上搜索。

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

2、.1691.613.903.651980122.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.

3、12529.9427.8620.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.76

4、20023656.573179.23410.8366.52200320042005200620072008200920102011201220134194.504892.185735.506706.197985.909905.1011484.1013606.8015988.3818331.3020796.503613.674333.185051.705898.797023.488890.3010312.8112207.1814320.87510.94496.10583.09678.83810.56826.10957.231147.991359.2769.8821.2250.9160.8969.

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

6、 0.0047 * z.diff.lag 0.02909 0.18067 0.161 0.8731Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1 Residual standard error: 0.08821 on 31 degrees of freedom Multiple R-squared: 0.2853, Adjusted R-squared: 0.2392 F-statistic: 6.188 on 2 and 31 DF, p-value: 0.005481 Value of test-statistic is: -3.046 Cri

7、tical values for test statistics:1pct 5pct 10pct tau1 -2.62 -1.95 -1.61从结果可以看出ADF统计量为-3.046,在1%, 5%, 10%的置信水平下均拒绝原 假设，认为序列x平稳。第二步模型的识别与定阶 模型的识别与定阶主要是通过对序列的自相关和偏自相关观察得到的。通过 R软件的到序列x的自相关和偏自相关图：图 4 序列 x 的自相关图图 5 序列 x 的偏自相关图通过自相关图和偏自相关图，可以发现序列 x 自相关明显拖尾，偏自相关有 结尾的性质，在第十阶的时候，PACF超出了两倍的标准差，这可能是由于偶然 因素引起的，可

8、以考虑使用AR(1 )和AR(ll)模型，但是为了数据的完整性， 尽量避免滞后期，可以考虑AR(1)模型。在 R 语言中有专门的 forecast 包，包里的 auto.arima 函数可以对序列模型进行智能识别建模。通过R语言的智能识别功能得到如下结果：Series: xARIMA(1,0,0) with zero meanCoefficients:ar10.4126s.e. 0.1512sigma人2 estimated as 0.006901: log likelihood=37.33AIC=-70.65AICc=-70.28BIC=-67.54通过如剑的输出结果我们发现结果与我们预测的

9、一样，均是AR(1)模型。 得到的模型方程为：X = 0.1426X “ + att 1t第三步模型的适应性检验模型的适应性检验实质上就是检验残差序列仇是否为白噪声序列。其中最主要 的是at序列的独立性检验。通过R软件得到模型的残差序列的时序图(偏)自相关图： 1- - v 、j-ResidualsOJ oCM oII51015202530下：图 6 残差序列的时序、自相关和偏自相关图通过对图形的观察，初步判断残差序列兔为白噪声。R 语言中的 Box.test 函数也可对残差序列进行白噪声检验。且得到的结果如m-EEdLAG LB p1,10.02128392 0.88400812,20.08

10、484764 0.95846353,30.43750962 0.93238944,40.696146050.95180335,52.02913305 0.84509946,62.03204959 0.91672757,72.24882031 0.94479578,82.45473056 0.96380859,98.50017079 0.484629310,10 10.925501510.363351511,11 10.939876320.448313912,12 13.844479050.3107497可以发现无论滞后几期， P 值均大于 0.05，所以接受原假设，认为残差序列 通过纯随机性检

11、验。故模型最终为ARIMA（l,0,0）,即为AR（1） 第四步模型参数的估计模型参数的估计主要有矩估计，最小二乘估计（ ML 估计），极大似然估，这里我们主要用的是最小二乘估计，估计的结果在第二步模型的识别中 R 语言智能建模中给出了，这里不做重复论述。模型的估计方程如下：存=01426X一 i + 兔第五步模型的预测将模型的拟合值与模型的真实值进行画图对比，得到如下图：1980 1985 1990 1995 2000 2005 2010Ti me图 7 拟合值真实值发现拟合值与真实值之间有一定差别，这是由于残差序列的存在，将拟合值与残差序列相加得到的结果与 X 序列作图，结果如下：图 8

12、加上残差的拟合值真实值观察图 8，可以发现拟合的特别好，所有的点都落在了真实值上。 我们用拟合的有效模型进行短期预测，比如我们预测2014年、2015年、2016、2017年和2018年的餐饮业零售总额。先预测2014年、2015年、2016 年、2017 年和2018年的x,再预测餐饮业零售总额。在R语言中可以使用predict函数对序列进行预测，这里我们进行向前5步预测, 得到点预测，区间预测。结果如下：$predTime Series:Start = 2014End = 2018Frequency = 11 -0.0254059516 -0.0104836538 -0.004326033

13、5 -0.0017851186 -0.0007366213 $seTime Series:Start = 2014End = 2018Frequency = 11 0.08307163 0.08986632 0.09097274 0.09115980 0.09119162 利用公式X Z x S.E.t 珞2logcaterts20l4=%oi4 + logcatertS20i3+ mean Cdifflogcaterts)可以得到预测值的区间估计，点估计和区间估计具体结果如下：序列X20142015201620172018点预测值-0.0254060-0.0104837-0.004326-0

14、.0017851-0.0007366下界值(L)-0.1882263-0.1866216-0.1826326-0.1804583-0.1794722上界值(U)0.13741440.16565430.17398050.17688810.1779990通过还原零均值、差分和对数变换，就可以得到江苏省餐饮业零售总额(caterts) 2014到2015年的预测值:利用公式：caterts20l4 = exp(logcaterts20l4)得到江苏省 2014 年餐饮业的零售手总额预测值，以此类推，分别得到2015,2016,2017,2018 年餐饮业零售总额的预测值。值得一体的是时间序列预测效果

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

16、.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=logged and diffed)x - scale(difflogcaterts,T,F) acf(x,10)#自相关pacf(x,1

17、0)#偏自相关 resid plotForecastErrors(fit$residuals)corr- acf(x,10) cov- acf(x,10,type = library(fUnitRoots) urdfTest(x) library(forecast) fit- auto.arima(x) fitresid- fit$residualcovariance)#ADF 单位根检验#建模#生成残差序列# 白噪声检验Result=0LAG=0LB=0p=0for(i in 1:12)Btest=Box.test(resid,type=Ljung-Box,lag=i) LAGi=iLBi=

18、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,lty=dashed) lines(L,col=blue,lty=dashed) plot(x,type=o)lines(fitted(fit),col=2,pch=7) point

19、s(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-

20、-0.0007366213+x2017+mean(difflogcaterts) forecast1 - exp(x2014) forecast2 - exp(x2015) forecast3 - exp(x2016) forecast4 - exp(x2017) forecast5 - exp(x2018)forecast1forecast2forecast3forecast4forecast5y - c(forecast1,forecast2,forecast3,forecast3,forecast4,forecast5) y - t(t(y)x_fit- fitted(fit)x_fit

21、- as.matrix(x_fit,ncol=1)i- 1z - 0for (i in 1:35) zi - x_fiti+logcatertsi+mean(difflogcaterts) x_fit_p- as.matrix(z,ncol=1)s - c(x2014,x2015,x2016,x2017,x2018)s - t(t(s)x_fit_p- rbind(x_fit_p,s)x_fit_pts- ts(x_fit_p,start=1979,frequency=1)#真实值与预测值图 plot.ts(caterts,main=(CATERING),ylab=CATERING,ylim=c(0,4500),xlim=c(197 8,2018),type=l,col=black)points(exp(x_fit_pts),pch=8,col=red)