時間序列ARIMA模型在R中的實現(xiàn)

1、。江M*fJIANGSUUNIVERSITY江蘇省餐飲業(yè)零售總額分析預(yù)測學(xué)校：江蘇大學(xué)學(xué)院：財經(jīng)學(xué)院班級：統(tǒng)計1201組員：韓亞瓊3120812015馬海燕3120812022顧君穎3120812020王培培3120812009陸金龍3120812029白卓3120812028完成時間：2014年12月13日星期六摘要引言數(shù)據(jù)分析原始數(shù)據(jù)的獲?。罕疚乃械臉颖緮?shù)據(jù)均來自江蘇統(tǒng)計年鑒一一2014»(得到的樣本數(shù)據(jù)參見表1:表1按行業(yè)分社會消費品零售總額年份社會消費品零售總額批發(fā)和零售業(yè)餐飲業(yè)其他行業(yè)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這里我們僅用到第三列數(shù)據(jù)，為了方便分析，我們將餐飲業(yè)零售總額序列命名為caterts第一步序列的平穩(wěn)性檢驗為判斷一個序列是否平穩(wěn)，我們主要通過時序圖以及自相關(guān)圖進行檢驗。對caterts做時序圖，有圖形發(fā)現(xiàn)有明顯的指數(shù)趨勢，序列非平穩(wěn)，也可以初步發(fā)現(xiàn)江蘇省的餐飲業(yè)零售總額逐年遞增，尤其是在新世紀以后，人們的生活水平逐年提高，對餐飲業(yè)的貢獻也增大：CATERING1980199020002010gooOTime圖1caterts序列時序圖因為原序列有明顯的指數(shù)趨勢，故先對數(shù)列進行對數(shù)變換得到新的

6、數(shù)列l(wèi)ogcatets,序列圖如下，具有明顯的非線性增長趨勢:圖2對數(shù)化后的時序圖對具有明顯線性趨勢的數(shù)列常用的平穩(wěn)化措施是差分，我們對logcaterts序列進行一階差分得到新的數(shù)列difflogcaterts,時序圖如下：loggedanddiffed圖3對數(shù)化和一階差分后的時序圖序列進行中心化處理，得到新的數(shù)列通過對時序圖分析發(fā)現(xiàn)數(shù)列具有平穩(wěn)性，為了方便分析，我們對difflogcatertsx。對x進行ADF檢驗（單位根檢驗）。R語言中有專門的fUnitRoots包，里面有urdftest功能，是專門進行序列的ADF單位根檢驗，通過R軟件得到如下結(jié)果：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從結(jié)果可以看出ADF

9、統(tǒng)計量為-3.046,在1%,5%,10%的置信水平下均拒絕原假設(shè)，認為序列x平穩(wěn)。第二步模型的識別與定階模型的識別與定階主要是通過對序列的自相關(guān)和偏自相關(guān)觀察得到的。通過R軟件的到序列x的自相關(guān)和偏自相關(guān)圖：圖4序列x的自相關(guān)圖Seriesx圖5序列x的偏自相關(guān)圖通過自相關(guān)圖和偏自相關(guān)圖，可以發(fā)現(xiàn)序列x自相關(guān)明顯拖尾，偏自相關(guān)有結(jié)尾的性質(zhì)，在第十階的時候，PACF超出了兩倍的標準差，這可能是由于偶然因素引起的，可以考慮使用AR(1)和AR(11)模型，但是為了數(shù)據(jù)的完整性，盡量避免滯后期，可以考慮AR(1)模型。在R語言中有專門的forecast包，包里的auto.arima函數(shù)可以對序列模

10、型進行智能識別建模。通過R語言的智能識別功能得到如下結(jié)果：Series:xARIMA(1,0,0)withzeromeanCoefficients:ar10.4126s.e.0.1512sigmaA2estimatedas0.006901:loglikelihood=37.33AIC=-70.65AICc=-70.28BIC=-67.54通過如劍的輸出結(jié)果我們發(fā)現(xiàn)結(jié)果與我們預(yù)測的一樣，均是AR(1)模型。得到的模型方程為：?=0.1426?一1+?第三步模型的適應(yīng)性檢驗?zāi)Ｐ偷倪m應(yīng)性檢驗實質(zhì)上就是檢驗殘差序列?是否為白噪聲序列。其中最主要的是?外序列的獨立性檢驗。通過R軟件得到模型的殘差序列的時

11、序圖(偏)自相關(guān)圖：Residuals圖6殘差序列的時序、自相關(guān)和偏自相關(guān)圖通過對圖形的觀察，初步判斷殘差序列?分為白噪聲。R語言中的Box.test函數(shù)也可對殘差序列進行白噪聲檢驗。且得到的結(jié)果如下: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可以發(fā)現(xiàn)無論滯后幾期，P值均大于0.05,所以接受原假設(shè)，認為殘差序列通過純隨機性檢驗。故模型最終為ARIMA(1,0,0),即為AR(1)第四步模型參數(shù)的估計模型參數(shù)的估計主要有矩估計，最小二乘估計（ML估計），極大似然估,這里我們主要用的是最小二乘估計，估計的結(jié)果在第二步模型的識別中R語言智能建模中給出了，這里不做重復(fù)論述。模型的估計方程如下：?=0.1426?一1+?第五步模型的預(yù)測將模型的擬合值與模型的真實值進行畫圖對比，得到如

13、下圖:19SC198519901995200020052010Time圖7擬合值一真實值發(fā)現(xiàn)擬合值與真實值之間有一定差別,這是由于殘差序列的存在，將擬合值與殘差序列相加得到的結(jié)果與X序列作圖，結(jié)果如下：圖8加上殘差的擬合值一真實值觀察圖8,可以發(fā)現(xiàn)擬合的特別好，所有的點都落在了真實值上。我們用擬合的有效模型進行短期預(yù)測，比如我們預(yù)測2014年、2015年、2016、2017年和2018年的餐飲業(yè)零售總額。先預(yù)測2014年、2015年、2016年、2017年和2018年的x,再預(yù)測餐飲業(yè)零售總額。在R語言中可以使用predict函數(shù)對序列進行預(yù)測，這里我們進行向前5步預(yù)測,得到點預(yù)測，區(qū)間預(yù)測。

14、結(jié)果如下：$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可以得到預(yù)測值的區(qū)間估計，點估計和區(qū)間估計具體結(jié)果如下:序列X20142015201620172018點預(yù)測值-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通過還原零均值、差分和對數(shù)變換，就可以得到江蘇省餐飲業(yè)零售總額(caterts(2014至U2015年的預(yù)測值：利用公式：?1?+?jrIoLjnIJ1IJ(bbk?1?exp7?0?得到江蘇省2014年餐飲業(yè)的零售手總額預(yù)測值，以此類推，分別得到2015,2016,2017,2018年餐飲業(yè)零售總額的預(yù)測值。值得一體的是時間序

16、列預(yù)測效果越好，則向前預(yù)測的階數(shù)就越少。所以在該模型中預(yù)測最好的是2014年。caterts20142015201620172018點預(yù)測值2088.2452474.9672951.4233528.5564222.97值得一體的是時間序列預(yù)測效果越好，則向前預(yù)測的階數(shù)就越少。所以在該模型中預(yù)測最好的是2014年。4、 參考文獻1王振龍.應(yīng)用時間序列分析M.北京：中國統(tǒng)計出版社，20102 PaulTeetor.R語言經(jīng)典實例M.北京：機械工業(yè)出版社，20133 RobertI.Kabacoff.R語言實戰(zhàn)M.北京：人民郵電出版社，20135、 附錄R語言代碼：'''rc

17、atering<-read.table("E:catering.txt",header=F)caterts<-ts(catering,start=1978,frequency=1)#數(shù)據(jù)時間序列化plot.ts(caterts,main=("CATERING"),ylab="CATERING")#序列圖logcaterts<-log(caterts)difflogcaterts<-diff(log(caterts)#對數(shù)差分plot(difflogcaterts,main="loggedanddiff

18、ed")x<-scale(difflogcaterts,T,F)#零均值化acf(x,10)#自相關(guān)pacf(x,10)#偏自相關(guān)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#接受原假設(shè)，認為resid序列為白噪聲# 預(yù)測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)# 預(yù)測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