时间序列分析课程报告

电子科技大学

时

间

序

列

分

析

基于时间序列分析的乘客人数预测

摘要：本文基于国际航空公司部分年限每个月的乘客数据，给出了一般非平稳时间序列的处理方式，通过对数据的预处理和分析，我们可以大致确定一个时间序列的平稳性，在确保了平稳性的基础之上后续的分析和建模才能顺利展开。通过构建ARIMA模型能够对大部分的时间序列产生较好的拟合结果，除此之外，本文还将数据集进行了一定的划分，使得训练和验证能够独立进行，能够检验模型的泛化能力，使得整个预测能够达到比较好的效果。

关键词：时间序列，ARIMA，分析和建模

一．研究背景

时间序列最早起源于古埃及人把尼罗河涨落的情况逐天记录下来,构成了所谓最早的时间序列相关的概念。初期的时间序列分析通常都是通过直观的数据比较或绘图观测,寻找序列中蕴含的发展规律，尼罗河周期性泛滥的规律发现就是依靠这种原始的分析方法。一个重要的转折点则是随机变量的引入，因为与时间序列相关的问题经常呈现出非常强的随机性，人们发现依靠单纯的描述性时序分析已不能准确地寻找出随机变量发展变化的规律。为了更准确地估计随机序列发展变化的规律，从20世纪20年代开始,一个大的趋势就是利用数理统计学的相关原理分析时间序列，研究的重心从表面现象的总结转移到分析序列值内在的相关关系上，由此开辟了应用统计学科[1]：时间序列分析。

关于时间序列分析方法，最早起源于1927年数学家Yule提出的自回归模型(AR 模型)，用来预测市场变化的规律。1931年，另一位数学家在AR模型的启发下,建立了移动平均模型(MA模型)，初步奠定了时间序列分析方法的基础。随后又逐渐提出其他模型，逐渐形成完整的时间序列分析方法。当前时间序列分析方法迈上了一个新的台阶,在工程领域方面的应用非常广泛，随着计算机技术和信号处理技术的迅速发展,时间序列分析理论和方法更趋完善。

具体而言，按照时间的顺序把随机事件变化发展的过程记录下来就构成了一个时间序列。对时间序列进行观察、研究，找寻它变化发展的规律，预测它将来的走势就是时间序列分析。简言之，时间序列分析的基本思想是根据系统有限长度的历史记录(即已知的观察数据)，建立能够比较精确地反映时间序列中所包含的动态依存关系的数学模型，并借以对系统的未来行为进行预测。典型的时间举例如下图所示：

图1.1951~1980美国劳工统计局统计的罢工人数

图2.1770~1869年间太阳黑子数统计

从模型构建的基础上讲，时间序列中的每个观察值大小，是由各种不同因素在同一时刻发生作用的综合结果。从这些影响因素发生作用的大小和方向变化的时间特性来看，这些因素造成的时间序列数据的变动分为下面几种：

趋势性：某个变量随着时间进展，呈现一种比较缓慢而长期的持续上升或者下降或者停留的不同属性的变动趋向，这种变动使得时间序列随机变量的方差在逐渐变化，进而凸显出一定的统计规律；

周期性：时间序列由于外部或者内部变量的影响随着时间的延伸，交替地出现高峰与低谷的规律；

随机性：由于随机白噪声的影响使得时间序列呈现出统计意义上的随机性；

因为时间序列本身的特点，常见的基本分析模型有下面几种：

1)P 阶自回归模型（AR(p)）:

即当前时刻的序列值由过去序列值的线性组合以及当前随机噪声构成；

2) Q 阶滑动平均模型（MA(q)）:

即当前时刻的序列值由过去的随机噪声和当前随机噪声构成；

3)ARMA模型(p,q):

即包含上述两种影响的时间序列问题，当前时刻值受过去的序列值以及过去的噪声和当前的噪声影响；

4)ARIMA模型(p,d,q):

由Box与Jenkins于上世纪七十年代提出，ARIMA的含义是单积自回归移动平均过程，其含义为：假设一个随机过程其经过d次差分后可以变换为一个平稳的自回归滑动平均过程，则该随机过程称为单积（整）自回归移动平均过程，其关键在于通过次数不多的差分使得非平稳的时间序列模型转化为平稳的时间序列模型(即普通ARMA模型)。

5)其他模型：

更为一般的非线性模型,如门限自回归、马尔科夫链；动态系统模型,如Kalman滤波，多层递阶预测；多元回归模型；变量场预测模型；

当然在实际应用中时间序列的分析方法远不止上述所列的，这里只是对基本的模型进行概述。时间序列在各个领域都有所应用，其中典型的应用领域就包括经济学的相关预测，信号处理领域，大数据挖掘以及气象信息预测等。举例而言，就经济学领域，诸如在证券领域中得到的观测数据列一般都具有较强的时间变化趋势，股票价格的数据都是以时间序列的形式出现的。因此,采用时间序列分析法对股市数据进行分析预测是可行的，并且很多文献都说明了它的有效性；除此之外间序列分析方法在数据挖掘中的应用也取得了一定的进展。研究者们利用数据挖掘对象，根据时间序列分析方法，提出了基于模糊集合的数据挖掘时间序列模式算法以及结合当下主流机器学习方法，比如支持向量机，随机森林等等；还有的则根据某些时间序列所具有的分形特征，分析了利用分形理论中的R/S分析，发现具有分形特

征的时间序列模式的方法；这些方法通过对大型数据库的海量数据分析提出了进行时间序列模式挖掘的算法，为用户的决策支持和趋势预测提供了依据。在气候预测领域，基于均生函数的时间序列预测延伸出了很多在气候分析及预测上有很好效果的模型，实验结果表明这类模型不但能较好地拟合历史实而且对未来1~5年的演变趋势也具有一定的预报能力[2]。

图3.时间序列分析在天气预测中的举例

二．基于R的乘客人数预测

本文希望通过借助于R平台，结合时间序列分析的基本方法，展现对于一般的非平稳模型的建模分析过程，借此加深对时间序列分析本身的理解以及相关理论知识的算法标准化，有助于后续对于时间序列分析的应用。

首先，本文的数据来源R本身的数据集“AirPassengers”,其中包含了1949-1960年每个月国际航空的乘客数量的数据，该数据集的特点为具有季节性和趋势性，是一种非平稳的时间序列，因此大致的思路是通过差分使得数据变得平稳，然后借助ARIMA模型来实现对规律的预测和分析。下图展示了该数据集：

图4.国际航空乘客数据集概览

首先，我对本文中涉及的数据进行时间序列的绘制，观察变量的总体趋势，并且拟合一条直线，通过这条拟合直线可以大致看出序列的趋势性；

图5.乘客人数从1949~1960的整体变化情况

从图上我们可以看出，随着时间轴的延伸，乘客人数是在逐渐递增的，这是大趋势，从拟合直线的走势可以看出来，这表明随着经济的逐步发展，越来越多的人们选择乘飞机去往世界各地，因此也符合现实情况的预期；除此之外，我们可以看到，在每一年的周期内部，乘客人数的变化也呈现出一定的规律性，大致是在每年的7，8月左右达到人数值的顶峰，这个我们在后面会做更精确的验证。总体而言，该数据满足具有季节性和趋势性的特点，是一种典型的非平稳数据。紧接着，我通过一个均值函数，求出这些年12个月份的乘客人数在每个月份点上的均值，并绘制如下的箱线图：

图7.乘客月平均人数箱线图

从上图可以看出，一年十二个月份中，每年的7，8月份是航空出行的旺季，因此我们也可以借此推测，每年的暑假阶段出行游览的人数较多，这也符合现实生活中的实际情况。因此，我们可以从简单的绘图以及对比，发现一些简单的推论，这些基于数据的合理推论也能给我们的现实生活产生一定的指导意义。

接下来我对这个时间序列进行相关的检验，首先我绘制出这个序列的自相关函数和偏相关函数的趋势图，如下图所示：

图8.未处理原始数据的自相关函数和偏相关函数

由自相关函数和偏相关函数，可以看出该时间序列不是平稳的，因为其自相关函数和偏相关函数并未呈典型的指数衰减，不能断定其是截尾的还是拖尾的。因此就上述分析而言，需要对原始数据进行差分，使得整体的序列呈现平稳性。为了使得数据更加符合正态性，在进行差分之前，我对原始的数据取对数之后在进行一阶差分；

图9.一阶差分之后的自相关函数和偏相关函数

在经过一阶差分之后，可以看出序列的衰减因而可以大致确定确定模型的阶数；但是因

为R内置了自动确定模型阶数的函数，因此我们这里直接调取函数。这里我根据序列的数据集大小采用了两种不同的处理方式，一种是选取包含原始所有数据的数据集，即通过原来所有的数据来预测下一个周期(即下一年)的乘客人数；另一种则是选取前十一年的数据作为样本集训练时间序列模型，剩下后一年的数据作为对比，相当于一个测试集来验证模型的预测效果。前者的处理是为了模型能够有充足的输入来拟合模型实现更精确的预测，而后一种则是借助训练集和测试集的划分来验证模型的预测性能，能够给出一定的泛化性能指标，这样即能够兼顾预测结果的完备性，也能够反映模型本身的性能优劣。

1）全序列：

原始样本包含了所有的十二年的数据，序列长度为144；

根据R自带的ARIMA模型参数自动匹配函数，最后有两组模型参数被选出来，即

分别对应时间序列分析模型的两组(p,d,q)，预测的置信区度为99.5%：1）第一个模型：(0,1,1)

其预测结果如下图所示：

模型参数检验如下及预测效果如下，蓝色部分为模型的预测下一个周期的乘客人数，：

图11.ARIMA(0,1,1)模型的下一个周期预测输出

2）第二个模型(2,1,1)

其模型参数检验结果和预测曲线如下所示：

图12. ARIMA(2,1,1) 模型的下一个周期预测输出

从预测曲线的直观看来，两个模型的预测差别较大，但是经过仔细的对比和挖掘，前者的信息量为AIC=1021，后者为AIC=1018.36，有所差别；除此之外，精确的预测值如下图显示：

图13.两种模型的预测值对比

我们可以看到，当阶数为(2,1,1)时，直观上来说更加符合7，8月份作为旺季的时候的规律，因此，虽然由R确定的最终模型有两个，但是明显地我们可以更倾向于选择ARIMA(2,1,1);整个模型的拟合度更高，而且从残差检验来说，ARIMA(2,1,1)的残差更加趋向于白噪声，说明其对趋势性和季节性的拟合更好，这点可以从下面的QQ-Plot可以看出，在QQ-Plot 中，残差拟合越接近理想的实对角线，说明模型拟合度越高，具体而言两种模型在数据中段的拟合程度都较好，只是在数据两端，ARIMA(2,1,1)比ARIMA(0,1,1)更加贴合，对比如下两图所示；

图14.ARIMA(2,1,1)的残差检验QQ-Plot

图15. ARIMA(0,1,1)的残差检验QQ-Plot

2）前11年的序列(除去最后一个周期的数据)：

这里是则是为了检验模型本身的性能，就如之前提到的，借助前11年的数据作为训练集，后一年的数据作为测试集的划分来验证模型的预测性能，能够给出一定的泛化性能指标，这样即能够兼顾预测结果的完备性，也能够反映模型本身的性能优劣。

值得一提的是，由于数据样本由原来的12个周期减少为11个，因此这里的ARIMA 需要重新确定阶数，显然地，当数据量改变的时候对于参数值的适应必然会有所改变。

由于减少了一个周期，因此，我先对这部分序列进行成分分解，观察其周期性和趋势性是否明显，结果如下图所示：

图16.前11个周期的序列分解图示

我们可以看到季节性和趋势性依然很强，说明序列整体的数学性质没有因为最后一年数据的剔除而发生明显改变。

因此根据上述分析，在构建模型之前，同样需要进行一阶差分处理；其自相关函数和偏相关函数如下所示：

图16.前11个周期的序列自相关函数和偏相关函数

根据自动定阶函数，确定此时的ARIMA模型为(1,1,0)，其残差检验如下图所示：

图17.前11个周期的残差检验QQ-Plot

由此可以看出模型的残差比较符合预期；对模型进行检验，结果如下：

接下来根据前11个周期建立的ARIMA模型来预测下一个周期的乘客人数，并且和实际的我们获取的数据真值进行比较，模型参数检验及结果如下图所示：

图18.前11个周期的ARIMA模型参数及其预测对比图

其中红色部分代表原始的真实数据，而蓝色部分代表由前11个周期预测的下一个周期的值，可以看到预测部分和其真实的分布重合度较好，说明ARIMA模型的拟合度较高，能够对原始数据集进行较好的拟合，从而印证了我们使用所有的数据集预测出的下一个周期乘客人数的可信度是较高的。

三．总结

本实验借助于R平台及其公开的数据集，对时间序列分析的建模过程进行了详细的探讨。对于一个给定的时间序列，我们在建立ARMA或者ARIMA模型时，需要做很多的预处理，这当中就包括：1)序列的平稳性判定；2)若是非平稳数据，则需要估计差分次数；3)自相关函数和偏相关函数的分析以及初步的定阶；4)模型残差的检验；做完这些步骤之后，我们就可以来建立对应的时间序列分析模型，这里对于不同的计算平台可能有不同的方式。我们可以根据自相关函数和偏相关函数给出一个比较粗略的阶次，当然也可以借助平台，比如R就提供了自动定阶的函数，可以通过自动定阶实现模型阶次的确定。通过后续的预测效果评估和分析，我们得到一个比较完备的模型，进而将其应用到实际的生产生活当中。

参考文献

[1]王燕1应用时间序列分析[M]1北京:中国人民大学出版社,2005:31

[2]张德宽,杨贤为,邹旭凯.均生函数)最优子集回归在高温极值预测中的应用[J].气,2003,29(4):44247.

代码

library("tseries")

library("forecast")

data(AirPassengers)

start(AirPassengers)

end(AirPassengers)

frequency(AirPassengers)

summary(AirPassengers)

boxplot(AirPassengers~cycle(AirPassengers),main="乘客人数的月均值对比")

opar<-par(mfrow=c(1,2))

acf(AirPassengers,main="自相关函数",lag.max=20)

pacf(AirPassengers,main="偏相关函数",lag.max=20)

par(opar)

dev.new()

y_cmp <- decompose(y)

plot(y_cmp)

airdf1 <- diff(AirPassengers,12)

Airdf1 <- diff(AirPassengers)

dev.new()

opar<-par(mfrow=c(1,2))

acf(airdf1,main="自相关函数",lag.max=20)

pacf(airdf1,main="偏相关函数",lag.max=20)

par(opar)

opar<-par(mfrow=c(1,2))

acf(Airdf1,main="自相关函数",lag.max=20)

pacf(Airdf1,main="偏相关函数",lag.max=20)

par(opar)

dev.new()

opar<-par(mfrow=c(1,2))

acf(b,main="自相关函数",lag.max=20)

pacf(b,main="偏相关函数",lag.max=20)

par(opar)

dev.new()

opar<-par(mfrow=c(1,2))

acf(y_df1,main="自相关函数",lag.max=20)

pacf(y_df1,main="偏相关函数",lag.max=20)

par(opar)

dev.new()

auto.arima(Airdf1,trace=T)#自动定阶

m1 <- arima(AirPassengers,order=c(2,1,1),seasonal=list(order=c(0,1,0),period=12)) dev.new()

qqnorm(m1$residuals)

qqline(m1$residuals)

Box.test(m1$residuals,type="Ljung-Box")#残差检验

dev.new()

y_frc1 <- forecast(m1,h=12,level=c(99.5))

plot(y_frc1)

lines(data,type="l",col="red")

plot(y_frc1,y_frc2)

lines(y_frc2,col="yellow")

adf.test(Airdf1,alternative="stationary",k=0)

auto.arima(b_df1,trace=T)

m1 <- arima(y,order=c(1,1,0),seasonal=list(order=c(0,1,0),period=12))

y_frc1 <- forecast(m1,h=12,level=c(99.5))

plot(y_frc1)

m011 <- arima(AirPassengers,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12)) dev.new()

y_frc2 <- forecast(m211,h=12,level=c(99.5))

plot(y_frc2)