铁路客流量预测
铁路客流量预测
目录
一、摘要........................................................................................................................................... 1 二、选题背景与意义 . ...................................................................................................................... 1 三、模型建立与求解 . ...................................................................................................................... 2
3.1、ARIMA 模型 .................................................................................................................. 2
3.1.1、自回归移动平均模型 . .......................................................................................... 2 3.1.2、季节性预测法 . ...................................................................................................... 3 3.1.3、模型求解 . .............................................................................................................. 4 3.2、灰色预测模型 . ................................................................................................................. 6
3.2.1、GM (1,1)模型 ...................................................................................................... 6
3.2.2、模型检验 . .............................................................................................................. 8
3.2.3、模型求解 . .............................................................................................................. 9
四、模型分析与结论 . .................................................................................................................... 11
4.1、方法分析 . ....................................................................................................................... 11 4.2、模型缺点 . ....................................................................................................................... 12 五、附录......................................................................................................................................... 12
一、 摘要
摘要：文章以铁路客流量的短期预测作为切入点，采用定量的时间序列分析方法，建立季节自回归综合移动平均（季节性ARIMA 模型）模型对时间序列进行量化分析。首先阐述基于该模型的预测的一般过程，即：平稳化处理、差分变换的阶数辨识、参数估计，时间序列模型的构建，然后利用标准BIC 值，确定较适合的季节自回归综合移动平均模型，取得了较为理想的预测效果。同时运用灰色预测模型建立铁路客流预测模型，对我国铁路客运量进行预测，灰色模型的方法简单，适合在数据少的情况下预测短期客流量，对未来的结果有很好的预测效果。
关键词：季节性ARIMA 灰色预测铁路客流量预测
二、 选题背景与意义
宏观上来讲铁路客流预测是铁路客运系统合理规划的基础，只有在对规划年度客流的流量、流向、流径进行合理预测与分析的基础之上，才能合理规划未来铁路客运系统的设施设备，合理安排运量，合理确定系统各阶段的发展目标使整个铁路客运系统与社会经济发展、生产力布局相适应，确保国民经济的正常发展。
微观层上来讲主要有以下三方面。
一是铁路客流量预测是铁路设备建设投资的重要依据。通过对各项客流预测结果分析，可以合理确定研究线路近期、中期、远期在路网中的功能和作用，从而为新线建设、旧线改造和相关客运场站技术设备修建与改造提供客观的依据。
二是铁路客流预测是编制铁路客流计划的基础。由于我国目前整体运能不足，再加上铁路运输自身的特点，在日常的客流运输组织中需要定期编制相应的客流计划，而准确的客流资料就是该项工作的基础，如果客流资料不完备就会造成运力资源分配的不平衡，从而致使客流滞塞及运力虚糜。
三是铁路客流预测是项目评价及投资估算的依据。铁路客运建设项目是否值得的投资，什么时候投资，投资规模如何，必须依据未来运量来确定。系统建成后，其寿命期内获利多少，也需要借助于逐年的未来运量才能估量和计算。如果没有科学、合理的运量为基础，就必然不能正确衡量和估算系统的经济成本和经济效益，致使经济评估失去真实性，导致投资决策的失误。
由以上分析可以看出铁路客流的预测对于系统的规划与建设、项目的投资与估算有着重要的依据。
三、 模型建立与求解
3.1、ARIMA 模型
随机时间序列分析模型可划分为3种不同类别：自回归模型(AR )[1]、滑动平均模型(MA )[2]以及自回归滑动平均模型(ARMA ) 。而自回归滑动平均模型研究的仅为平稳时间序列，而对于非平稳时间序列则通常采用自回归综合移动平均模型ARIMA 。ARIMA 模型亦可分为带趋势性的模型ARIMA (p , d , q )，和既带有趋势又有季节性趋势的模型
ARIMA (p , d , q )⨯(P , D , Q ) s 。
自回归移动平均过程是由自回归和移动平均两部分组成的随机过程，形式化表示为
ARMA (p , q ) ，其中p 和q 分别为自回归和移动平均部分的最大阶数。ARMA (p , q ) 的数学
表达式为：
X t =φ1X t -1+φ2X p -2+ +φp X t -p +εt +θ1εt -1+θ1εt -2+ +θq εt -q
提取公因式，得到如下式子:
(1-φ1L -φ2L 2- -φp L p ) X t =(1+θ1L +θ2L 2+ +θq L q ) εt
将其中的乘积项替换，亦可表示为:
Φ(L ) X t =Θ(L ) εt
其中，Φ(L ) 和Θ(L ) 分别表示自变量L 的p ，q 阶特征多项式。
3.1.1、 自回归移动平均模型
ARMA 即自回归综合移动平均模型，它满足如下条件，{x t }为自回归整和移动平均
序列，记为ARIMA (p , d , q )，其中，d 为整和阶数，p 为自回归系数，q 为移动平均系数。
在一般的自回归移动平均模型中，无季节性，仅有趋势性。假设x t 表示随机序列，并假定：
Lx t =x t -1
其中L 是滞后算子。如果存在非负整数d ，满足：
φ(L )∇d x t =θ(L )εt
式中函数表示为：
⎧Φ(L )=1-Φ1L -Φ2L 2- -Φp L p ⎪⎪2p ⎨Θ(L )=1+Θ1L +Θ2L + +Θp L ⎪d d ∇=1-L ()⎪⎩
φ(L ) 与θ(L ) 互质，且L ≤1，存在E (εt ) =0，{εt }是白噪音序列，φp θp ≠0；E (εt 2)
3.1.2、 季节性预测法[3]
某些不平稳的时间序列既具有趋势演化性，又会随进行周期性的演化，通常若一个序
列的演化周期为S ，那么该序列将每隔S 个时间间隔均呈类似的变化。假定有整数D ≥0，以及随机序列{x t , t =0, ±1,... }，满足式：
s
φ(L s )∇d s x t =θ(L )εt
则时间序列{x t }表示季节性ARIMA (P , D , Q )过程，其中∇s =1-L S ，为季节差分算子，S 为季节性周期，则：
⎧∇s x t =(1-L S )x t =x t -x t -s ⎪
⎨D
S D D -1D -1
⎪⎩∇s x t =(1-L )x t =x t -x t -s
其中，D 为季节性差分阶数。且：
⎧Φ(L s )=1-Φ1L s -Φ2L 2s - -Φp L ps ⎪⎨s s 2s ps ΘL =1+ΘL +ΘL + +ΘL ()⎪12p ⎩
其中P 为季节性自回归阶数，(Φ1, Φ2, , Φn )为季节性自回归部分的参数，Q 为季节性滑动平均阶数，Θ1, Θ2, , Θp 为季节性移动平均阶数部分的参数。将两式融合，变为一般的季节ARIMA 模型，即：
()
φp (L )ΦP (L s )∇D ∇d x t =θq (L )ΦQ (L s )εt
这里，p 、d 、q 、P 、D 和Q 的不同是为了调整不同算子的阶数，可称得到的季节ARIMA
模型为ARIMA (p , d , q )⨯(P , D , Q ) 。
s
3.1.3、 模型求解
我们从国家统计局得到的2008.1-2016.9铁路客流量月数据作为时间序列数据，用上述模型进行分析，并通过建立的模型来预测未来一年铁路客流量的变化情况。
从上面的时间序列图可以看出，在每年快春节的时候和十一的时候，客流量是明显高于其他每月的，这也正与实际相吻合，受到春节，假日的影响；同时从图中可以很直观的看出整个客流量呈现出稳定的上升趋势，所以说铁路客流量具有明显的周期性和趋势性，所以我们采用季节性ARIMA 模型，即求出p 、d 、q 、P 、D 和Q 的值则确定了模型。 由于时间序列明显有上升趋势，所以该时间序列是非平稳的，所以我们先进行一阶差分处理，消除其显著的趋势性，得到下图。
从一阶差分序列图可以发现序列图围绕0值上下波动，其方差明显有界，
所以时间序列
的趋势性有所消除，而一阶处理后的铁路客流量自相关和偏自相关函数值如下所示。
如图所示，ACF 与PACF 均呈拖尾形态在零值邻域波动，而且1，2，10，12阶相关函数大于0，与春节，国庆等假日很有关。为了取得更好的效果，使时间序列更加合理，我们再比对二阶非季节性差分处理的结果，以求得更恰当的参数。
由二阶差分序列图可以看出效果并没有很大的改善，在2012年12月之前的序列是更加平稳了，但后面时间的并不理想，所以我们还是先采用一阶差分处理，即选取d =1，从图取得拖尾阶数选择p =2，q =2。
下图为一阶季节性差分和一阶非季节性差分的自相关图和偏自相关图。
由于在实际情况中[]，p , d , q , P , D , Q ∈(0,1,2)且不全为0，所以也验证了上面选取p , q
也是合理的。由于一般情况下，季节性差分阶数D =1，由于季节自回归阶数P ，季节移动平均阶数Q 难以确定，为精确起见，我们同时建立多个模型，在系数显著的情况下使用了 BIC 准则来进行比较。我们考虑对p , d , q , P , D , Q 取不同的值共有9种组合，来算BIC 与考察序列残差是否是白噪声。在这9种不同的组合中我们选取BIC 的值最小的组合。下面是我们得到的表。
p 2222
d 1111
q 2222
P 0001
D 1111
Q 0120平稳的R 方 0.721 0.725 0.799 0.722
标准化BIC 15.011 15.31 15.054 15.32
22222
11111
22222
11222
11111
12012
0.741 0.797 0.767 0.710 0.808
15.309 15.125 15.205 15.485 15.134
由该表，我们得到了p =2, d =1, q =2, P =0, D =1, Q =0的组合，此时BIC=15.011最小。因此我们选用参数定阶对客流量进行预测，经SPSS 处理后得到未来一年铁路客流量的变化以及与原数据比较得到的残差的自相关和偏自相关图。
由预测时序图可以看出整个趋势以及每月的变化预测的还是较为合理。
残差序列的样本自相关函数与偏自相关函数基本可控制数均可控制95%的置信区间之内，因此，残差序列为白噪声过程（随机变化过程）。在季节性ARIMA 预测法在短期内能输出较理想的预测结果，但随预测时间的增加，预测的误差将逐渐增大，因为预测时间的增加使得预测置信区间的宽度也变大，所以该模型更适用于短期预测。
3.2、灰色预测模型
灰色系统预测理论的基本思路是按某种规则将已知的数据序列构成非动态的或动态的白色模块，然后按照某种变换解决来求解未来的灰色模型。在灰色系统理论中，常用的模型是微分方程所描述的动态方程，最简单的是基于灰色系统理论模型GM
(1,1)模型的预测分
析。灰色预测分析可分为几类，即数列预测，灾变预测，季节性灾变预测，拓扑预测及系统综合预测。
GM (1,1)模型[4]
灰色理论的微分方程模型称为GM 模型，GM (1,1)表示一阶、单个变量的微分方程。
GM (1,1)是一阶单序列的线性动态模型，用于时间t 序列预测的是其离散形式的微分方程模
型，具体形式为
dx
+ax =u dt
由上式可知，这是一个单变量x 对时间的一阶微分方程，是连续的，实际使用的是其离散的单个数据形式。 设有数列x
(0)
共有n 个观察值x (0)(1)，x (0)(2)，x (0)(3)， ，x (0)(n ) ，对x
(1)
(0)
作一次累加
生成，得到新的数列x ，表达式为
x (i ) =∑x (1)(m ), i =1,2, , n
(1)
m -1
i
対一阶生成数列x 建立预测模型，其方程为
(1)
dx (1)
+ax (1)=u dt
式中：a ，u 为待估参数，分别称为发展灰数和内生控制灰数。 将上式的离散形式展开，可得
k =1, x (1)(2)=
⎛1⎫ a -(x (1)(1)+x (1)(2))⎪+u ;
⎝2⎭k =2, x (1)(3)=
⎛1⎫ a -(x (1)(2)+x (2)(3))⎪+u ;
⎝2⎭ k =n , x (1)(n ) =
⎛1⎫ a -(x (1)(n -1) +x (2)(n )) ⎪+u ;
⎝2⎭
将两个待估模型参数表示为向量形式得
⎡a ⎤ˆ=⎢⎥ a
⎣u ⎦
将上述离散方程组用最小二乘法求解，得
ˆ=(B T B ) -1B T y n a
ˆ代入上式，解微分方程，得到GM (1,1)的预测模型为 将a
u ⎫u ⎛
ˆ(1)(k +1) = x (0)(1)-⎪e -ak + x
a ⎭a ⎝
式中
y n =[x (0)(2),x (0)(3) , x (0)(n )]T ;
⎡1(1)⎤(1)
-(x (1)+x (2)) 1⎢2⎥⎢⎥⎢-1(x (1)(2)+x (1)(3)) 1⎥
⎥ B =⎢2
⎢ ⎥⎢⎥⎢1(1)⎥(1)-(x (n -1) +x (n )) 1⎢⎥⎣2⎦
模型检验
灰色预测模型的检验，有关联检验、后验检验和残差检验。残差检验分两种：一是相对
误差，二是绝对误差。检验步骤为设原始序列：
x 0=(x (0)(1),x (0)(2), x (0)(n ))
灰色预测模型序列为：
ˆ0=(x ˆ(0)(1), x ˆ(0)(2) , x ˆ(0)(n ) ) x
计算残差
ˆ(0)(n ) ε(0)(n ) =x (0)(n ) -x
计算相对误差
∆=
ε(0)(n )
x (0)(n )
计算x (0)(n ) 的均值和方差为
1n (0)1n (0)2
=∑x (k ), S 1=∑[x (k ) -]2
n k =1n k =1
计算ε(0)(n ) 的均值和方差为
1n (0)1n (0)2
=∑ε(k ), S 2=∑[ε(k ) -]2
n k =1n k =1
称C =
S 20
为均方差比值也叫后验差比，称p =P (ε(0)-
率。指标C 越小越好，p 越大越好。一般地，将模型精度等级分为四级，如下表：
模型精度等级 1级(好) 2级(合格) 3级(勉强) 4级(不合格)
C C
C ≥0.65
p >0.95 0.80
p ≤0.70
如果关联度、方差、小误差概率和相关误差比都在允许范围之内时，则可用所建模型进行预测，否则应进行残差修正。
模型求解
用2008～2016年的数据来预测2016年11和12月及2017年1～10月各月客运量，建立铁路客流量灰色预测模型。
首先将每年各月的数据提取出来，将11月的数据提取出来为 2008 2009 2010 2011 2012
客运量/万人 13215 13480 14761 15828 16600
将2008～2015年各年11月的客流量能够得到
年份/年
2013 17972 2014 19471 2015 21231
x (0)=[x (0)(1),x (0)(2),x (0)(3),x (0)(4),x (0)(5),x (0)(6),x (0)(7),x (0)(8)]=
[13215，13480,14761,15828,16600,17972,19471,21231]
求得一次累加生成数列
x (1)=[13215,26695,41456,57284,73884,91856,111327,132558]
经MATLANB 处理后求得
⎡-19955 1⎤⎡26695⎤⎢-34075.5 1⎥⎢41456⎥⎢⎥⎢⎥⎢-49370 1⎥⎢57284⎥⎢⎥⎢⎥B =⎢-65584 1⎥，y n =⎢73884⎥
⎢-82870 1⎥⎢91856⎥⎢⎥⎢⎥-101591.5 1111327⎢⎥⎢⎥⎢-121942.5 1⎥⎢132558⎥⎣⎦⎣⎦
⎡-0.0735⎤ˆ=⎢ a ⎥⎣ 12054⎦
即
a =-0.0735, u =12054
所以
x (0)(1)=13215,
于是可以得到预测模型为
u
=-163910 a
ˆ(1)(k +1) =177125e 0.0735k -163910(k =0,1,2 ) x
以表格形式列出预测值和实际值
k
1
13215
2
26724
3
41263
4
56912
5
73753
6
91880
7
111390
8
132390
9
154980
ˆ(1)(k ) x
x (0)(k )
实际值
13215 13215
13509 13480
14539 14761
15648 15828
16842 16600
18126 17972
19509 19471
20997 21231
22598
计算绝对误差序列和相对误差序列分别为
∆(0)={0,29,222,180,242,154,38,234}
Φ={0,0.22%,1.5%,1.13%,1.46%,0.86%,0.2%,1.1%}
由程序运行后得到p=1，C=0.24
由同样的方法，可以得到2016.10～2017.9一年内的数据，列在下表中：
时间 预测值 时间 预测值
2016.10 22378 2017.4 25617
2016.11 22598 2017.5 24858
2016.12 25374 2017.6 25643
2017.1 22380 2017.7 27060
2017.2 23272 2017.8 28041
2017.3 24355 2017.9 26282
将前几年的数据和下一年度的数据统计绘制在下图中
由图中可以看出，数据预测效果还是符合客流量整体的趋势，且还是有相应的峰值出现，这与客流量周期性，季节性相符合。但是相比12月之前的数据，预测出来的数值偏大，将12月单独罗列拿出来看
可以看出2015年12月的客流量是比较低的，但是对整个预测趋势来说还是呈现出上升态势的，可以说灰色预测对整体的把握还行，但一旦出现小幅波动之类的情况，预测结果的可信度就不是很强。
四、 模型分析与结论
下面是ARIMA 预测法和灰色预测法对接下来一年所作预测的时序图，可以看出两种方法的预测效果还是比较接近，只是对个别的值灰色预测法还有所欠缺，下面是两种方法的分析总结。
4.1、方法分析
1. 由ARIMA 模型得到的拟和结果可知短期时间序列的预测精度是比较高的。由此可见，
自回归时间序列预测法是一种重要的预测方法，其模型比较简单，对资料的要求比较单一，只需变量本身的历史数据，在实际中有着广泛的适用性。在应用中，应根据所须解决的问题及问题的特性等因素来综合考量并选择相对优化的模型。 2. 灰色预测方法简单，虽然该模型是建立在高等数学基础上，但计算步骤简单，可以借助
计算机软件很容易计算出来，计算时间短。 3. 灰色预测模型需要的数据少。由于灰色预测把随机过程看作灰色过程，所以预测只根据
实际情况选择适量的数据即可。 4. 灰色预测可有效的处理贫信息和数据少的情况。在一定时间段内预测的精度较高，但是
随着信息的增加，不断进入灰色系统时，会发现预测效果越来越差，灰色系统不适合长期的预测，不能用该模型预测未来的所有值。
4.2、模型缺点
1. 在选择参数时有很强的主观因素，从自相关图和偏自相关图中确定p , d , q , P , D , Q ，的值还缺乏一定的科学性。并且未进行更多次的实验选取，从中选择出更优的方案。 2. 在灰色预测过程中，没有对异常值处理做过多的分析说明，对12月数据的变化还不能
有效的解释和说明。
五、 附录
文献
[1]张志雷. 自相关过程的ARIMA 控制图[J].统计与决策，2012，(6)
[2]张立杰，寇纪淞，李敏强等. 基于自回归移动平均及支持向量机的中国棉花价格预测[J].统计与决策，2013，(6)
[3]肖良，基于季节性ARIMA 模型的居民消费水平预测[A].安徽，宿州，2016 [4]黄召杰，冯硕. 灰色预测模型在铁路客流预测中的应用[A].交通科技与经济，第16卷1期，2014.02
代码
function GM(x0) %灰色系统GM （1，1）预测 %x0=[12794.1483800000 13804.1483800000 15508.1483800000 16480.1483800000 15310.1483800000 16631.1483800000 18143.1483800000 22910.1483800000]; %10月份 %x0=[12794 13804 15508 16480 15310 16631 18143 22910];36.91 %11月数据 %x0=[11639.8671300000 12199.8671300000 13495.8671300000 14452.8671300000 16121.8671300000 22365.8671300000 23733.8671300000 19506.8671300000]; %12月份 %x0=[11907.6900400000 13289.6900400000 12731.6900400000 15202.6900400000 16475.6900400000 18764.6900400000 19057.6900400000 17857.6900400000 21168.6900400000];%1月份 %x0=[13243.1692100000 13984.1692100000 14613.1692100000 16115.1692100000
15966.1692100000 14437.1692100000 16368.1692100000 19683.1692100000 14682.8046300000 22124.8046300000 15338.0531400000 20884.0531400000 15616.5531400000 21526.5531400000 15342.0948000000 20880.0948000000 15877.8723300000 22493.8723300000 15216.0129600000 22893.0129600000 15781.4348300000 21445.4348300000
24505.1692100000];%2月份 %x0=[12425.8046300000 12370.8046300000 14660.8046300000 15027.8046300000 17424.8046300000 18624.8046300000 21812.8046300000];%3月份 %x0=[11415.0531400000 12283.0531400000 13062.0531400000 16245.0531400000 17295.0531400000 19636.0531400000 23693.0531400000];%4月份 %x0=[11970.5531400000 13195.5531400000 14091.5531400000 15184.5531400000 16539.5531400000 19344.5531400000 23193.5531400000];%5月份 %x0=[11732.0948000000 11785.0948000000 13630.0948000000 16492.0948000000 18309.0948000000 19722.0948000000 23466.0948000000];%6月份 %x0=[11511.8723300000 11905.8723300000 13718.8723300000 15701.8723300000 17648.8723300000 20103.8723300000 24535.8723300000];%7月份 %x0=[11413.0129600000 12361.0129600000 13554.0129600000 15871.0129600000 17641.0129600000 20869.0129600000 25361.0129600000];%8月份 %x0=[12139.4348300000 11822.4348300000 13462.4348300000 16557.4348300000 18840.4348300000 20629.4348300000
23561.4348300000]; %9月份
x0=[36.55 45.46 46.33 45.13 46.33 43.5 44.88 44.17 43.9 43.58 43 42.38 43.1 42.33 42.4 43.5... 42.88 41.48 41.68 42.72 41.63 41.93 42.7 42.4 42.28 42.23 41.22 42.73 42.02 42.25 42.3 42.65];
%x0=[13215 13480 14761 15828 16600 17972 19471 21231];
NUM=32;
%[31.26 32.09 33.39 35.45 40.52 43.52 47.86 51.45 55.6 60.14 64.82 68.65 73.22 80.22 87.69 93.97 99.16 103.38 109.46 114.6 119.85 124.92 132.04 139.45 150.93 165.13 181.93 196.83 212.38 228.07 241.45 251.03];
T=input('请输入T ：'); %预测接下来几年 x1=zeros(1,length(x0)); B=zeros(length(x0)-1,2); yn=zeros(length(x0)-1,1); hatx0=zeros(1,length(x0)+T); hatx00=zeros(1,length(x0)); hatx1=zeros(1,length(x0)+T); epsilon=zeros(length(x0),1); omega=zeros(length(x0),1); for i=1:length(x0) for j=1:i
x1(i)=x1(i)+x0(j);%累加生成 end end x1
for i=1:length(x0)-1
B(i,1)=(-1/2)*(x1(i)+x1(i+1)); B(i,2)=1; yn(i)=x0(i+1); end
hatA=(inv(B'*B))*B'*yn ; %求a ，u for k=1:length(x0)+T
hatx1(k)=(x0(1)-hatA(2)/hatA(1))*exp(-hatA(1)*(k-1))+hatA(2)/hatA(1); end
hatx0(1)=hatx1(1); for k=2:length(x0)+T
hatx0(k)=hatx1(k)-hatx1(k-1);%累减还原 end
hatx0 %历史数据的模拟值 for i=1:length(x0)%开始检验 epsilon(i)=x0(i)-hatx0(i);
omega(i)=(epsilon(i)/x0(i))*100; end
test1=hatA(2)/hatA(1) %hatA(1)为a ，hatA(2)为u test2=hatA(1)
c=std(epsilon)/std(x0); p=0;
for i=1:length(x0)
if abs(epsilon(i)-mean(epsilon))
p=p/length(x0) if p>0.95&c
disp('预测精度好')
elseif p>0.85&c
disp('预测精度合格')
elseif p>0.70&c
disp('预测精度勉强合格')
elseif p0.65 disp('预测精度不合格') end
for i=1:length(x0)
hatx00(i)=hatx0(i); end
z=1:length(x0)+T;y=1:length(x0);
plot(y+2007,x0,'r-',z+2007,hatx0,'b-','linewidth',2); hold on;
%scatter(y+1984,x0,'r'); %scatter(z+1984,hatx0,'b'); %axis([1975 2020 0 500]);
%text(1990,440,'红色：实际值','color','red') %text(1990,410,'蓝色：模拟值','color','blue') grid
for k=1:NUM
u(k)=abs(hatx0(k)-x0(k)) end
v=sum(u(1:NUM))/NUM