一、作業要求
自選時間序列完成時間序列的建模過程，要求序列的長度>=100。
報告要求以下幾部分內容：
- 數據的描述：數據來源、期間、數據的定義、數據長度。
- 作時間序列圖並進行簡單評價。
- 進行時間序列的平穩性檢驗，得出結論，不平穩時間序列要進行轉化，最終平穩。
- 進行自相關、偏自相關圖，得出模型的階數。
- 對時間序列模型進行擬合，得出參數的估計值。
- 檢驗模型的殘差項，判斷模型是否合格，給出模型最終的估計結果。
- 應用建立的時間序列模型進行預測。
二、數據描述
數據來源：國家統計局——統計數據——月度數據——交通運輸——旅客運輸量
時間范圍選擇“2005-”，表示2005年至今
點擊“下載”，格式選擇CSV，並重命名為“旅客運輸量.csv”
http://data.stats.gov.cn/easyquery.htm?cn=A01
本次使用的數據為表中的第8行——鐵路客運量當期值(萬人)
期間：2005年1月至2019年10月
其中2005年至2017年的數據用來建立模型
2018年和2019年的數據用於和預測結果比較
數據的定義：數據類型為時間序列（ts）
#載入必要的R程序包
library(fUnitRoots)
## Warning: package 'fUnitRoots' was built under R version 3.5.3
## Loading required package: timeDate
## Warning: package 'timeDate' was built under R version 3.5.3
## Loading required package: timeSeries
## Warning: package 'timeSeries' was built under R version 3.5.3
## Loading required package: fBasics
## Warning: package 'fBasics' was built under R version 3.5.3
library(forecast)
## Warning: package 'forecast' was built under R version 3.5.3
#讀入旅客運輸量.csv"
transport<-read.csv("旅客運輸量.csv",header=F)
#提取鐵路客運量當期值(萬人)月度數據(2005年1月至2019年10月)
tr<-transport[8,180:3]
tr<-t(tr)
tr<-as.numeric(tr)
#轉換成時序數據
tr<-ts(tr,frequency=12,start=c(2005,1))
此處可以查看全部數據，如下所示。
tr#查看數據 ## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec ## 2005 9300 10600 9300 9100 9700 8600 10800 11200 9400 10000 8600 8500 ## 2006 10700 11300 9900 9900 10700 9600 12000 12200 10200 11000 9273 9200 ## 2007 9900 11090 11973 10255 11400 10200 13100 13495 11448 12079 10300 10700 ## 2008 11900 12850 11855 11622 11663 11466 13794 14059 12496 12570 10800 10333 ## 2009 13282 13591 11800 12490 12888 11519 14188 15007 12179 13580 11065 10893 ## 2010 12724 14220 14090 13269 13784 13364 16001 16200 13819 15284 12346 12189 ## 2011 15195 15722 14112 15545 15309 15076 18160 17862 16138 16256 13413 13146 ## 2012 16468 15573 14457 16452 14877 16226 17984 18517 16914 15086 14185 14815 ## 2013 18757 14044 16854 17502 16232 18043 19931 20287 19197 16407 15557 17375 ## 2014 19050 15975 18054 19843 19037 19456 22386 23515 20986 17919 17056 22427 ## 2015 17850 19290 21554 21091 21219 20614 24776 25539 21802 22686 18816 18200 ## 2016 21161 24112 21242 23900 22886 23200 26818 28007 23918 25001 20409 20768 ## 2017 24756 25525 22624 26504 26397 24077 29378 30692 24884 27621 22703 23219 ## 2018 24564 26081 27612 28900 26827 27834 32276 34340 28253 30467 25177 25164 ## 2019 28342 29112 27860 30536 30801 30735 35570 37884 29873 31903
數據長度：178
n<-length(tr) #此處預留2018年和2019年的數據進行驗證 tr_pred<-window(tr,start=c(2018,1),end=c(2019, 10)) tr<-window(tr,start=c(2005,1),end=c(2017, 12)) tr_len<-length(tr)
三、時間序列圖
#畫出時序圖 plot.ts(tr) abline(lm(tr~time(tr)))
圖1 原始時間序列圖
如上圖1所示，可以看到鐵路客運量逐年上升，明顯是非平穩序列，呈拋物線式增長趨勢，后面我們將利用二次函數對其進行擬合。又可以觀察到在每年之中有明顯的周期性，后面我們將提取這一周期性趨勢並剔除。
通過簡單移動平均可以更好地看出這一趨勢，如下圖2和圖3所示。
#通過簡單移動平均進行平滑處理並觀察 plot.ts(ma(tr,3),main="Simple Moving Averages (k=3)")
圖2 移動平均（k=3）
plot.ts(ma(tr,5),main="Simple Moving Averages (k=5)")
圖3 移動平均（k=5）
四、平穩性檢驗
利用ADF檢驗判斷原始序列的平穩性。
#檢驗平穩性 adfTest(tr) ## ## Title: ## Augmented Dickey-Fuller Test ## ## Test Results: ## PARAMETER: ## Lag Order: 1 ## STATISTIC: ## Dickey-Fuller: 0.08 ## P VALUE: ## 0.6398 ## ## Description: ## Mon Feb 17 08:07:03 2020 by user: lenovo
可以看到，P VALUE大於0.05，原始序列不平穩。
因前面觀察得知該序列具有周期性，故下面畫出月度圖進行進一步觀察，如下圖4和圖5所示。
#繪出月度圖 monthplot(tr)
圖4 月度圖（1）
seasonplot(tr,year.labels="TRUE")
圖5 月度圖（2）
至此，以年為周期的趨勢更加明顯，7月、8月鐵路客運量較多，11月和12月則較少。
為了分離周期性趨勢，我們進行季節性分解，如下圖6所示。
為了采用乘法形式分解，我們首先對原始數據取了對數。
#進行季節性分解 ltr<-log(tr) fit<-stl(ltr,s.window="period") plot(fit)
圖6 季節性分解
fit$time.series#展現分解結果 ## seasonal trend remainder ## Jan 2005 0.0029836818 9.146746 -0.0119599765 ## Feb 2005 0.0157985766 9.150486 0.1023242527 ## Mar 2005 -0.0246972016 9.154227 0.0082399529 ## Apr 2005 0.0008667277 9.158901 -0.0437380821 ## May 2005 -0.0017750290 9.163575 0.0180810274 ## Jun 2005 -0.0383998121 9.168771 -0.0708539884 ## Jul 2005 0.1357229189 9.173967 -0.0223889054 ## Aug 2005 0.1586769877 9.179102 -0.0141094555 ## Sep 2005 0.0079087099 9.184236 -0.0436793923 ## Oct 2005 0.0143873586 9.191987 0.0039658463 ## Nov 2005 -0.1457697038 9.199739 0.0055485025 ## Dec 2005 -0.1257032021 9.208776 -0.0352516340 ## Jan 2006 0.0029836818 9.217814 0.0572014651 ## Feb 2006 0.0157985766 9.225411 0.0913487074 ## ......
利用分解結果，即可剔除周期性因素，如下圖7所示。
#剔除周期性因素 m<-exp(fit$time.series[,1]) tr_m<-tr/m #繪圖 plot.ts(tr_m)
圖7 剔除周期性后時間序列
此時趨勢性仍十分明顯，下面采用二次函數擬合並剔除趨勢性，如下圖8所示。
#二次函數回歸 x1<-1:length(tr) x2<-x1^2 lm<-lm(tr_m~x1+x2) summary(lm) ## ## Call: ## lm(formula = tr_m ~ x1 + x2) ## ## Residuals: ## Min 1Q Median 3Q Max ## -3447.6 -448.3 49.3 531.7 5106.0 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 9.760e+03 2.295e+02 42.529 < 2e-16 *** ## x1 2.589e+01 6.749e+00 3.836 0.000182 *** ## x2 5.179e-01 4.164e-02 12.438 < 2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 943.2 on 153 degrees of freedom ## Multiple R-squared: 0.9652, Adjusted R-squared: 0.9647 ## F-statistic: 2121 on 2 and 153 DF, p-value: < 2.2e-16 #剔除趨勢項 tr_m2<-tr_m-(lm$coefficients[1]+lm$coefficients[2]*x1+lm$coefficients[3]*x2) plot.ts(tr_m2)
可以看到，相關系數為0.9652，各因素均顯著（***）。
回歸表達式為：\[y=9.76\times {{10}^{3}}+25.89x+0.5179{{x}^{2}}\]
圖8 剔除周期性和趨勢項后時間序列
對此時剔除周期性和趨勢項后的序列進行平穩性檢驗和白噪聲檢驗。
#平穩性檢驗 adfTest(tr_m2) ## Warning in adfTest(tr_m2): p-value smaller than printed p-value ## ## Title: ## Augmented Dickey-Fuller Test ## ## Test Results: ## PARAMETER: ## Lag Order: 1 ## STATISTIC: ## Dickey-Fuller: -11.8908 ## P VALUE: ## 0.01 ## ## Description: ## Mon Feb 17 08:07:03 2020 by user: lenovo #白噪聲檢驗 Box.test(tr_m2, type="Ljung-Box") ## ## Box-Ljung test ## ## data: tr_m2 ## X-squared = 4.755, df = 1, p-value = 0.02921
P VALUE、p-value均小於0.05，序列平穩，且非白噪聲。
五、自相關、偏自相關圖
下面繪出自相關和偏自相關圖，如圖9和圖10所示。
#繪出自相關和偏自相關圖 acf(tr_m2,lag.max=48)
圖9 自相關圖
如果樣本(偏)自相關系數在最初的d階明顯大於2倍標准差范圍，而后幾乎95％的(偏)自相關系數都落在2倍標准差的范圍以內，而且由非零自相關系數衰減為在零附近小值波動的過程非常突然。這時通常視為(偏)自相關系數截尾，截尾階數為d。
ACF拖尾。若認為4階截尾，其后還有4個數值超出二倍標准差，$\frac{4}{48-4}$=9%>5%，不可。
pacf(tr_m2,lag.max=48)
圖10 偏自相關圖
PACF為3階截尾。
在4階及以后的45個數值中，有2個超出2倍標准差。
$\frac{2}{45}$=4.4%<5%，符合要求。
故選用AR(3)模型。
六、模型擬合、參數估計
下面利用arima函數進行ARIMA（3,0,0）擬合。
#進行ARIMA模型擬合 tr_m2<-c(tr_m2) arima_tr_m2<-arima(x=tr_m2,order=c(3,0,0),method="ML") arima_tr_m2 ## ## Call: ## arima(x = tr_m2, order = c(3, 0, 0), method = "ML") ## ## Coefficients: ## ar1 ar2 ar3 intercept ## -0.1522 -0.1707 0.2552 1.9318 ## s.e. 0.0772 0.0767 0.0771 64.8146 ## ## sigma^2 estimated as 749431: log likelihood = -1276.63, aic = 2563.27
擬合結果為${{y}_{t}}^{\prime }={{y}_{t}}-1.9318$，
${{y}_{t}}^{\prime }=-0.1522{{y}_{t-1}}^{\prime }-0.1707{{y}_{t-2}}^{\prime }+0.2552{{y}_{t-3}}^{\prime }$。
下面利用該表達式進預測。
#給出下一個預測值 a=forecast(arima_tr_m2,h=1) a$mean[1] ## [1] 310.1413 #手動計算 arima_tr_m2$coef[1]*(tr_m2[tr_len]-arima_tr_m2$coef[4])+ arima_tr_m2$coef[2]*(tr_m2[tr_len-1]-arima_tr_m2$coef[4])+ arima_tr_m2$coef[3]*(tr_m2[tr_len-2]-arima_tr_m2$coef[4])+ arima_tr_m2$coef[4] ## ar1 ## 310.1413
自動計算值和手動計算值一致，可以驗證結果的正確性。
利用此預測值，可以計算2018年1月鐵路客運量的預測值。
趨勢項$y=9.76\times {{10}^{3}}+25.89\times 157+0.5179\times {{157}^{2}}=26590$，
周期項為1.0029881，
故${{\widehat{y}}_{2018-1}}=(310.1413+26590)\times 1.0030=26981.0$。
七、殘差檢驗
下面對殘差進行白噪聲檢驗和正態分布檢驗，繪出QQ圖，如下圖11所示。
error<-arima_tr_m2$residuals #檢驗殘差是否為白噪聲 Box.test(error, type="Ljung-Box") ## ## Box-Ljung test ## ## data: error ## X-squared = 0.2025, df = 1, p-value = 0.6527
p-value大於0.05，殘差是白噪聲。
#繪制殘差QQ圖 qqnorm(error) qqline(error)
圖11 殘差QQ圖
通過QQ圖可以看到，大部分殘差數據落在直線上，符合正態分布，但存在少數極端值不符合。
#殘差正態性檢驗 shapiro.test(error) ## ## Shapiro-Wilk normality test ## ## data: error ## W = 0.93947, p-value = 3.23e-06
通過Shapiro-檢驗，p-value小於0.05，可知殘差不滿足正態分布。
八、結果預測
采用循環方式逐一預測，生成2018-1至2019-10的預測序列，並與原序列進行對比，如下圖12所示。
#逐一重復預測
preds<-NULL
newdata<-NULL
for(i in 1:22)
{
ts0=c(tr_m2,preds)
arima_tr_m2<-arima(x=ts0,order=c(3,0,0),method="ML")
a=forecast(arima_tr_m2,h=1)
preds=c(preds,a$mean[1])
}
#增加趨勢項
x1<-(tr_len+1):(tr_len+22)
x2<-x1^2
preds2<-preds+(lm$coefficients[1]+lm$coefficients[2]*x1+lm$coefficients[3]*x2)
#增加周期項
mm<-rep(m,length.out=22)
preds3<-preds2*mm
#繪圖展示
plot(c(tr,tr_pred),type='l',col='darkgreen',lwd=2)
lines((tr_len+1):(tr_len+22),preds3,col='red')
points((tr_len+1):(tr_len+22),preds3,col='red',pch=20)
圖12 預測結果
上圖中綠色線為原始數據，紅色點線為預測數據，二者數值相近，趨勢一致，預測效果較好。
下面利用平均絕對百分比誤差來評估預測效果。
#計算平均絕對百分比誤差(mean absolute percentage error) (MAPE<-mean(abs(preds3-tr_pred)/tr_pred)) ## [1] 0.03639169
可以看到，平均絕對百分比誤差約為3.6%，在可接受范圍內。
九、自動預測
利用R函數auto.arima全自動進行建模和預測，效果如下圖13所示。
#自動選擇模型預測
preds_auto<-NULL
for(i in 1:22)
{
ts0=c(tr,preds_auto)
ts0<-ts(ts0,frequency=12,start=c(2005,1))
arima_tr_auto<-auto.arima(ts0)
a=forecast(arima_tr_auto,h=1)
preds_auto=c(preds_auto,a$mean[1])
}
arima_tr_auto
## Series: ts0
## ARIMA(2,1,1)(0,1,2)[12]
##
## Coefficients:
## ar1 ar2 ma1 sma1 sma2
## -0.4743 -0.4567 -0.6502 -0.2436 -0.1732
## s.e. 0.0860 0.0798 0.0800 0.0768 0.0736
##
## sigma^2 estimated as 827511: log likelihood=-1349.4
## AIC=2710.8 AICc=2711.33 BIC=2729.39
plot(c(tr,tr_pred),type='l',col='darkgreen',lwd=2)
lines((tr_len+1):(tr_len+22),preds_auto,col='red')
points((tr_len+1):(tr_len+22),preds_auto,col='red',pch=20)
圖13 自動預測結果
(MAPE_auto<-mean(abs(preds_auto-tr_pred)/tr_pred)) ## [1] 0.03936244
可以看到，R為我們選擇了具有周期項的(2,1,1)(0,1,2)[12]模型，並自動進行了預測，平均絕對百分比誤差約為3.9%。此時誤差稍大，可能是由於進行了兩次差分（包含一次季節差分）而損失了部分信息量。
十、收獲感悟
通過本次大作業，我對時間序列分析有了更深的了解和認識，對ARIMA模型有了更清晰的掌握，並通過編程實現熟悉了R中關於時間序列的命令，以及對時間序列建模和預測的過程，收獲頗豐。
十一、參考資料
- 時間序列分析與應用. Jonathan D.Cryer, Kung-Sik Chan著，潘紅宇等譯，機械工業出版社
- 應用時間序列分析. 白曉東著，清華大學出版社
- R語言預測實戰. 游皓麟著，電子工業出版社
- R語言實戰. Robert I.Kabacoff著，王小寧，劉擷芯，黃俊文等譯，人民郵電出版社
- R數據科學. 哈德利·威克姆，加勒特·格羅勒芒德著，陳光欣譯，人民郵電出版社
- https://blog.csdn.net/zxy_clover/article/details/79750267
- https://blog.csdn.net/tantaixf/article/details/83148901
- https://blog.csdn.net/c1z2w3456789/article/details/80752541
- https://blog.csdn.net/desilting/article/details/39013825