時間序列作業(yè)

1、 時間序列分析基于R期末考察報告班 級 1314132專 業(yè) 數(shù)學與應用數(shù)學學 號 131413213姓 名 盧巖巖 指導教師 劉常勝 時 間 2016.06.201、 問題重述一個塑料制造商生產(chǎn)的產(chǎn)品A的月銷售量數(shù)據(jù)，如下表所示。1 742 741 896 951 1030 697 700 793 861 1032 776 774 885 93815 1126 898 932 1055 1109 1285 1030 1099 1204 1274 1468 1107 1223 132629 1422 1637 1165 1290 1303 1486 1611 1216 1349 1436 15

2、55 1608 1208 134143 1473 1604 1528 1131 1296 1453 1600 1420 971 1066 1170 1403 1119 78357 901 1023 1209 1013繪制該序列時序圖，自相關圖，首先直觀考察該序列的特點。使用因素分解方法擬合該序列的發(fā)展，并預測下一月該產(chǎn)品的銷售量。然后對數(shù)據(jù)進行差分運算后，使用ARIMA模型對觀察序列建模，并進行預測。 二、函數(shù)表示read.csv(" ", header= ) 讀入數(shù)據(jù)ts(x,start= , freq= ) 指定為時間序列變量plot() 繪制時序圖acf(x) 繪制序

3、列自相關圖pacf(x) 繪制序列偏自相關圖for(i in 1: 2) print(Box.test(x, lag=6*i ) ) 白噪聲檢驗forecast(x, h= ,level= ) 預測 decompose(x ,type= ) 確定性因素分解HoltWinters(x ,alphy = ,bata= ,gamma= ,seasonal= ) HoltWinters指數(shù)平滑模型diff( ) 差分運算arima(x, order=c(p,d,q), method=” ”) 擬合ARIMA模型 三擬合并預測3.1獲得觀察值序列：（1）將數(shù)據(jù)文件讀入R并存入數(shù)據(jù)框x.讀入命令和結(jié)果如下

4、>x<-read.csv("G:/8.csv",head=T) xMonth Mduct.A.for.a.plastics.manufacturer1 1月1日 7422 1月2日 7413 1月3日 8964 1月4日 9515 1月5日 10306 1月6日 6977 1月7日 7008 1月8日 7939 1月9日 86110 1月10日 103211 1月11日 77612 1月12日 77413 2月1日 88514 2月2日 93815 2月3日 112616 2月4日 89817 2月5日 93218 2月6日

5、 105519 2月7日 110920 2月8日 128521 2月9日 103022 2月10日 109923 2月11日 120424 2月12日 127425 3月1日 146826 3月2日 110727 3月3日 122328 3月4日 132629 3月5日 142230 3月6日 163731 3月7日 116532 3月8日 129033 3月9日 130334 3月10日 148635 3月11日 161136 3月12日 121637 4月1日 134938 4月2日 143639 4月3日 155540 4月4日 160841 4月5日 120842 4月6日 13414

6、3 4月7日 147344 4月8日 160445 4月9日 152846 4月10日 113147 4月11日 129648 4月12日 145349 5月1日 160050 5月2日 142051 5月3日 97152 5月4日 106653 5月5日 117054 5月6日 140355 5月7日 111956 5月8日 78357 5月9日 90158 5月10日 102359 5月11日 120960 5月12日 1013（2）提取第二列數(shù)據(jù)> x1<-x,2x1 1 742 741 896 951 1030 697 700 793 861 1032 776 774 88

7、5 93815 1126 898 932 1055 1109 1285 1030 1099 1204 1274 1468 1107 1223 132629 1422 1637 1165 1290 1303 1486 1611 1216 1349 1436 1555 1608 1208 134143 1473 1604 1528 1131 1296 1453 1600 1420 971 1066 1170 1403 1119 78357 901 1023 1209 1013（3）將以上數(shù)據(jù)轉(zhuǎn)化為時間序列數(shù)據(jù)x2<-ts(x1,freq=12,start=1)3.2平穩(wěn)性檢驗畫出時序圖和自相

8、關圖>plot(x2,type="o",main="一個塑料制造商的產(chǎn)品A的月銷售量",xlab="年份",ylab="月銷售量") > acf(x2)注：其中圖中的年份表示實際數(shù)據(jù)的月份數(shù)據(jù)，12個月是實際每月的12天。這樣數(shù)據(jù)代表的意義與實際數(shù)據(jù)相吻合。（以下所有年份代表的數(shù)據(jù)都是月份數(shù)據(jù)）結(jié)果如下圖所示:月度銷售量序列時序圖月銷售量序列自相關圖時序圖呈現(xiàn)出逐年變化的趨勢，及以年為周期規(guī)則的周期性，并且自相關圖衰減向零的速度較慢，長期位于零軸的一邊，該序列一定不是平穩(wěn)序列。由于該非平穩(wěn)序列有顯著的趨

9、勢及固定的周期變化，可進行確定性因素分解，本文采用乘法模型：式中，代表序列的長期趨勢波動；代表序列的季節(jié)性（周期性）變化；代表隨機波動。3.3確定性因素分解x2.fit<-decompose(x2,type="mult")x2.fit$x Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec1 742 741 896 951 1030 697 700 793 861 1032 776 7742 885 938 1126 898 932 1055 1109 1285 1030 1099 1204 12743 1468 1107

10、1223 1326 1422 1637 1165 1290 1303 1486 1611 12164 1349 1436 1555 1608 1208 1341 1473 1604 1528 1131 1296 14535 1600 1420 971 1066 1170 1403 1119 783 901 1023 1209 1013$seasonal Jan Feb Mar Apr May Jun Jul1 1.0729463 0.9930701 0.9931157 0.9826870 0.9562760 1.0928037 0.93477772 1.0729463 0.9930701 0.

11、9931157 0.9826870 0.9562760 1.0928037 0.93477773 1.0729463 0.9930701 0.9931157 0.9826870 0.9562760 1.0928037 0.93477774 1.0729463 0.9930701 0.9931157 0.9826870 0.9562760 1.0928037 0.93477775 1.0729463 0.9930701 0.9931157 0.9826870 0.9562760 1.0928037 0.9347777 Aug Sep Oct Nov Dec1 1.0331357 0.978215

12、2 1.0000889 0.9992614 0.96362222 1.0331357 0.9782152 1.0000889 0.9992614 0.96362223 1.0331357 0.9782152 1.0000889 0.9992614 0.96362224 1.0331357 0.9782152 1.0000889 0.9992614 0.96362225 1.0331357 0.9782152 1.0000889 0.9992614 0.9636222$trend Jan Feb Mar Apr May Jun Jul1 NA NA NA NA NA NA 838.70832 9

13、14.5417 952.0833 979.6250 989.4583 1010.0833 1048.7500 1093.87503 1267.6667 1270.2083 1281.7917 1309.2917 1342.3750 1356.9167 1349.54174 1393.5000 1419.4167 1441.8750 1436.4583 1408.5417 1405.2917 1425.62505 1328.1667 1279.2083 1218.8750 1188.2500 1180.1250 1158.1667 NA Aug Sep Oct Nov Dec1 852.8750

14、 870.6667 878.0417 871.7500 882.58332 1125.2083 1136.2917 1158.1667 1196.4167 1241.08333 1358.2917 1385.8333 1411.4167 1414.2500 1393.00004 1435.4167 1410.4167 1363.5000 1339.3333 1340.33335 NA NA NA NA NA$random Jan Feb Mar Apr May Jun Jul1 NA NA NA NA NA NA 0.89285042 0.9019071 0.9920829 1.1573873

15、 0.9235568 0.9648848 0.9205308 1.08456483 1.0793021 0.8775922 0.9607473 1.0306042 1.1077519 1.1039600 0.92348804 0.9022502 1.0187430 1.0859328 1.1391418 0.8968380 0.8732129 1.10532275 1.1227665 1.1178078 0.8021586 0.9129231 1.0367513 1.1085223 NA Aug Sep Oct Nov Dec1 0.8999750 1.0109201 1.1752383 0.

16、8908214 0.91007762 1.1053831 0.9266442 0.9488291 1.0070822 1.06527493 0.9192620 0.9611673 1.0527493 1.1399616 0.90589044 1.0816058 1.1074943 0.8294092 0.9683608 1.12498315 NA NA NA NA NA$figure 1 1.0729463 0.9930701 0.9931157 0.9826870 0.9562760 1.0928037 0.9347777 8 1.0331357 0.9782152 1.0000889 0.

17、9992614 0.9636222$type1 "multiplicative"attr(,"class")1 "decomposed.ts"> plot(x2.fit)結(jié)果如下圖所示：因素分解擬合綜合輸出圖示3.4進行HoltWinters三參數(shù)指數(shù)平滑并繪制擬合效果圖> b<-HoltWinters(x2)bHolt-Winters exponential smoothing with trend and additive seasonal component.Call:HoltWinters(x = x2)Sm

18、oothing parameters: alpha: 0.07424144 beta : 0.5194062 gamma: 0.1735239Coefficients: ,1a 970.571884b -42.331143s1 38.713698s2 -15.936653s3 34.078847s4 -63.295181s5 -80.085910s6 55.488371s7 -94.214048s8 -62.216597s9 -47.603478s10 33.357650s11 0.978706s12 -53.359675> plot(b) HoltWinters三參數(shù)指數(shù)平滑擬合效果圖

19、3.5做1年期預測，并繪制預測效果圖> library(zoo)8> library(forecast)> x2.forec<-forecast(b,h=12)> x2.forec Point Forecast Lo 80 Hi 80 Lo 95 Hi 95Jan 6 966.9544 710.31598 1223.5929 574.45979 1359.4491Feb 6 869.9729 611.70685 1128.2390 474.98904 1264.9569Mar 6 877.6573 616.48613 1138.8285 478.23046 127

20、7.0841Apr 6 737.9521 472.27151 1003.6327 331.62869 1144.2756May 6 678.8303 406.75549 950.9050 262.72781 1094.9327Jun 6 772.0734 491.49948 1052.6473 342.97262 1201.1742Jul 6 580.0398 288.70979 871.3699 134.48897 1025.5907Aug 6 569.7061 265.28049 874.1318 104.12728 1035.2850Sep 6 541.9881 222.10855 86

21、1.8677 52.77453 1031.2017Oct 6 580.6181 242.96005 918.2762 64.21467 1097.0215Nov 6 505.9080 148.21972 863.5963 -41.12902 1052.9451Dec 6 409.2385 29.36657 789.1104 -171.72548 990.2025上表的數(shù)據(jù)就是下一月前12天該制造商的產(chǎn)品銷售量的預測值。第一列是預測時期；第二列就是進行點預測每天產(chǎn)品銷售量的數(shù)值； plot(x2.forec)月度銷售量HoltWinters三參數(shù)指數(shù)平滑序列預測圖確定性因素分解方法可以對觀察值序

22、列很好的解釋，但不能充分提取觀察值序列的有效信息，導致模型擬合精度通常不夠精確。為了彌補確定性因素分解方法的不足，可使用ARIMA模型對觀察值序列建模。首先做1階12步差分，并繪制差分后序列時序圖plot(x2.fit)> x3<-diff(diff(x2),12)> plot(x3)產(chǎn)品月銷售量差分后序列時序圖繪制差分后序列自相關圖和偏自相關圖> acf(x3)> pacf(x3)產(chǎn)品月銷售量差分后序列自相關圖產(chǎn)品月銷售量差分后序列偏自相關圖產(chǎn)品A月銷售量的原始時序圖顯示序列既有趨勢也有周期特征。1階12步差分后序列時序圖呈現(xiàn)出平穩(wěn)特征。差分后序列自相關圖和偏自

23、相關圖都顯示出拖尾屬性，所以首先嘗試擬合加法季節(jié)模型ARIMA（1，（1,12），1）模型。擬合ARIMA（1，（1,12），1）模型> c1<-arima(a,order=c(1,1,1),seasonal=list(order=c(0,1,0),period=12)> c1Call:arima(x = x2, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 0), period = 12)Coefficients: ar1 ma1 -0.4826 1.0000s.e. 0.1327 0.1188sigma2 esti

24、mated as 68825: log likelihood = -329.89, aic = 665.79殘差序列白噪聲檢驗> for(i in 1: 2) print(Box.test(c1$residual,lag=6*i) Box-Pierce testdata: c1$residualX-squared = 58.676, df = 6, p-value = 8.357e-11 Box-Pierce testdata: c1$residualX-squared = 112.27, df = 12, p-value < 2.2e-16通過上面殘差序列的白噪聲檢驗顯示，殘差序

25、列非白噪聲，該模型擬合效果不理想。這說明加法季節(jié)模型并不適合擬合這個序列。所以，需要考慮乘積季節(jié)模型。綜合前面的差分信息，擬合的乘積模型為。擬合模型c2<-arima(x2,order=c(8,1,1),seasonal=list(order=c(1,1,0),period=12)> c2Call:arima(x = x2, order = c(8, 1, 1), seasonal = list(order = c(1, 1, 0), period = 12)Coefficients: ar1 ar2 ar3 ar4 ar5 ar6 ar7 ar8 -0.1129 0.0685 -

26、0.2174 -0.0119 0.9275 0.0072 -0.1007 0.1582s.e. 0.6635 0.1847 0.1599 0.0727 0.0545 0.6127 0.1579 0.1765 ma1 sar1 -0.0268 -0.5743s.e. 0.6679 0.1476sigma2 estimated as 3357: log likelihood = -266.47, aic = 554.94殘差序列白噪聲檢驗for(i in 1: 2) print(Box.test(c2$residual,lag=6*i) Box-Pierce testdata: c2$residualX-squared = 2.6794, df = 6, p-value = 0.8479 Box-Pierce testdata: c2$residualX-squared = 5.331, df = 12, p-value = 0.946殘差序列的白噪聲檢驗中，