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

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

6、數(shù)列l(wèi)ogcatets,序列圖如下，具有明顯的非線性增長(zhǎng)趨勢(shì):圖2對(duì)數(shù)化后的時(shí)序圖對(duì)具有明顯線性趨勢(shì)的數(shù)列常用的平穩(wěn)化措施是差分，我們對(duì)logcaterts序列進(jìn)行一階差分得到新的數(shù)列difflogcaterts,時(shí)序圖如下：loggedanddiffed圖3對(duì)數(shù)化和一階差分后的時(shí)序圖序列進(jìn)行中心化處理，得到新的數(shù)列通過對(duì)時(shí)序圖分析發(fā)現(xiàn)數(shù)列具有平穩(wěn)性，為了方便分析，我們對(duì)difflogcatertsx。對(duì)x進(jìn)行ADF檢驗(yàn)（單位根檢驗(yàn)）。R語(yǔ)言中有專門的fUnitRoots包，里面有urdftest功能，是專門進(jìn)行序列的ADF單位根檢驗(yàn)，通過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)計(jì)量為-3.046,在1%,5%,10%的置信水平下均拒絕原假設(shè)，認(rèn)為序列x平穩(wěn)。第二步模型的識(shí)別與定階模型的識(shí)別與定階主要是通過對(duì)序列的自相關(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ì)，在第十階的時(shí)候，PACF超出了兩倍的標(biāo)準(zhǔn)差，這可能是由于偶然因素引起的，可以考慮使用AR(1)和AR(11)模型，但是為了數(shù)據(jù)的完整性，盡量避免滯后期，可以考慮AR(1)模型。在R語(yǔ)言中有專門的forecast包，包里的auto.arima函數(shù)可以對(duì)序列模

10、型進(jìn)行智能識(shí)別建模。通過R語(yǔ)言的智能識(shí)別功能得到如下結(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ù)測(cè)的一樣，均是AR(1)模型。得到的模型方程為：?=0.1426?一1+?第三步模型的適應(yīng)性檢驗(yàn)?zāi)Ｐ偷倪m應(yīng)性檢驗(yàn)實(shí)質(zhì)上就是檢驗(yàn)殘差序列?是否為白噪聲序列。其中最主要的是?外序列的獨(dú)立性檢驗(yàn)。通過R軟件得到模型的殘差序列的時(shí)

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

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

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

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

17、atering<-read.table("E:catering.txt",header=F)caterts<-ts(catering,start=1978,frequency=1)#數(shù)據(jù)時(shí)間序列化plot.ts(caterts,main=("CATERING"),ylab="CATERING")#序列圖logcaterts<-log(caterts)difflogcaterts<-diff(log(caterts)#對(duì)數(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單位根檢驗(yàn)urdfTest(x)library(forecast)#建模fit<-auto.arima(x)fitresid<-fit$residual#生成殘差序列residplotForecastErrors(fit$residuals)#白噪聲檢驗(yàn)Result=0

19、LAG=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è)，認(rèn)為resid序列為白噪聲#預(yù)測(cè)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",l

20、ty="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ù)測(cè)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