一、作业要求
自选时间序列完成时间序列的建模过程，要求序列的长度>=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
.nav-tabs {
display: inline-table;
max-height: 500px;
min-height: 44px;
overflow-y: auto;
background: white;
border: 1px solid #ddd;
border-radius: 4px;
}
.tabset-dropdown > .nav-tabs > li.active:before {
content: "";
font-family: 'Glyphicons Halflings';
display: inline-block;
padding: 10px;
border-right: 1px solid #ddd;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open > li.active:before {
content: "";
border: none;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open:before {
content: "";
font-family: 'Glyphicons Halflings';
display: inline-block;
padding: 10px;
border-right: 1px solid #ddd;
}
.tabset-dropdown > .nav-tabs > li.active {
display: block;
}
.tabset-dropdown > .nav-tabs > li > a,
.tabset-dropdown > .nav-tabs > li > a:focus,
.tabset-dropdown > .nav-tabs > li > a:hover {
border: none;
display: inline-block;
border-radius: 4px;
}
.tabset-dropdown > .nav-tabs.nav-tabs-open > li {
display: block;
float: none;
}
.tabset-dropdown > .nav-tabs > li {
display: none;
}
-->