时间序列

2018-09-03

2018-09-03
时间序列

1 介绍

  时间序列分析是数理统计这一学科的一个重要分支。时间序列分析的最基本的理论基础是20世纪40年代分别由Norbort Wiener和Andrei Kolmogonov独立给出的.在60年代后期，时间序列分析在Box-Jenkins提出一套比较完善的建模理论及方法之后便迅速发展起来．近30年来，计算技术和计算机的发展又赋予时间序列分析以新的活力，使之成为自然科学、社会科学领域中不可缺少的数据分析工具。

  通俗地讲，按照时间的顺序把随机事件变化发展的过程记录下来就构成了一个时间序列；对时间序列进行观察、研究，找寻它变化发展的规律，预测它将来的走势就是时间序列分析。

1.1 时间序列分析的目的

  (1)预测序列的未来发展情况。

  (2)分析序列的基本趋势、季节和随机项的组成。

  (3)分析特定的数据集合、拟合理论模型，尤其是建立数学模型，进而进行模型结构分析和实证研究

1.2 时间序列分析的主要应用

  (1)预测

  (2)政策因素分析

  (3)信号的提取

  (4)随机系统描述

2 定义

  严格来讲，时间序列是指将同一统计指标的数值按其发生的时间先后顺序排列而成的数列。表现形式为\(X_1,X_2,...X_t\),记作\(\left\{X_t\right\}\)。拿到一个序列观察值之后，首先要对它的平稳性和纯随机性进行检验，这两个重要的检验称为序列的预处理。根据检验的结果可以将序列分为不同的类型，对不同类型的序列我们会采用不同的分析方法。

3 平稳序列

3.1 平稳序列介绍

  首先我们来看看平稳时间序列，平稳时间序列有两种定义，根据限制条件的严格程度，分为严平稳时间序列和宽平稳时间序列。

  • 严平稳是一种条件比较苛刻的平稳性定义，它认为只有当序列所有的统计性质都不会随着时间的推移而发生变化时，该序列才能被认为平稳。我们知道随机变量族的统计性质完全由它们的联合概率分布族决定。所以严平稳的核心就是分布函数在参数平移中的不变性。

  • 宽平稳是使用序列的特征统计量来定义的一种平稳性，它认为序列的统计性质主要由它的低阶矩决定，所以只要保证序列低阶矩平稳（二阶），就能保证序列的主要性质近似稳定。宽平稳的核心是矩在参数平移中的不变性。

  那么，严平稳和宽平稳的联系和区别是什么呢？可以看出，严平稳是对序列联合分布的要求，以保证序列所有的统计特征都相同；而宽平稳只要求序列二阶平稳，对于高于二阶的矩没有任何要求。所以通常情况下，严平稳序列也满足宽平稳条件，而宽平稳序列不能反推严平稳成立。

  但这不是绝对的。比如服从柯西分布的严平稳序列就不是宽平稳序列，因为它不存在一、二阶矩，所以无法验证它二阶平稳，严格地讲，只有存在二阶矩的严平稳序列才能保证它一定也是宽平稳序列。宽平稳一般推不出严平稳，但当序列服从多元正态分布时，二阶平稳可以推出严平稳。

  在实际应用中，研究最多的是宽平稳随机序列，以后见到平稳随机序列，如果不是特别注明，指的都是宽平稳随机序列，如果随机序列不满足平稳条件，就称为非平稳序列。

  判别平稳序列的三个标准：

  (1)均值是与时间t无关的常数。

  (2)方差是与时间t无关的常数。

  (3)协方差是只与时期间隔k有关，与时间t无关的函数。

以上三个标准同时符合，则为平稳序列。

  我们为什么要研究平稳时间序列呢？数理统计学的基本原理，是利用样本信息来推测总体信息，根据数理统计学常识，显然要分析的随机变量越少越好，而每个变量获得的样本信息越多越好。因为随机变量越少，分析的过程就会越简单，而样本容量越大，分析的结果就会越可靠。

  但是时间序列分析的数据结构有它的特殊性。对随机序列\(\left\{...,X_1,X_2,...,X_t,...\right\}\)而言，它在任意时刻\(t\)的序列值\(X_t\)都是一个随机变量，而且由于时间的不可重复性，该变量在任意一个时刻只能获得唯一的样本观察值。由于样本信息太少，如果没有其他的辅助信息，通常这种数据结构是没有办法进行分析的。而序列平稳性概念的提出可以有效地解决这个问题。

  在平稳序列场合，序列的均值等于常数，意味着原本含有可列多个随机变量的均值序列\(\left\{\mu_t,t\in T\right\}\)变成了一个常数序列\(\left\{\mu,t\in T\right\}\)，原本每个随机变量的均值\(\mu(t\in T)\)只能依靠唯一的一个样本观察值\(x_i\)去估计\(\hat\mu=x_i\)，现在由于\(\mu_t=\mu(\forall t\in T)\)，于是每一个样本观察值\(x_t(\forall t\in T)\)，都变成了常数均值\(\mu\)的样本观察值\[\hat\mu=\bar X=\frac{\sum_{i=1}^nx_i}{n}\]  这极大地减少了随机变量的个数，并增加了待估参数的样本容量。换句话说，这大大降低了时序分析的难度，同时也提高了对均值函数的估计精度。同理，根据平稳序列二阶矩平稳的性质，可以得到基于全体观察样本计算出来的延迟k自协方差函数估计值、总体方差估计值和延迟\(k\)自相关系数估计值。

3.2 平稳序列建模

  拿到一个序列之后，首先要判断它的平稳性，常用的检验方法是根据时序图和自相关图显示的特征做出判断的图检验方法。通过平稳性检验，序列可以分为平稳序列和非平稳序列两大类。对于平稳序列，我们有一套非常成熟的平稳序列建模方法。

  但是，并不是所有的平稳序列都值得建模。只有那些序列值之间具有密切相关关系，历史数据对未来发展有一定影响的序列，才值得我们花时间去挖掘历史数据中的有效信息，用来预测序列未来的发展。如果序列值之间没有任何相关性，那就意味着该序列是一个没有记忆的序列，过去的行为对将来的发展没有丝毫影响，这种序列叫做纯随机序列（又叫做白噪声序列），纯随机序列没有任何分析价值。

  ARMA模型的全称是自回归移动平均模型，它是目前最常用的拟合平稳序列的模型。它又可以细分为AR模型（自回归模型）、MA模型（移动平均模型）和ARMA模型。

3.2.1 AR模型

  对于偏自相关系数p阶截尾，自相关系数拖尾大的数据，我们采用AR模型。设\(\left\{X_t\right\}\)为零均值平稳序列，若满足如下的\(p\)阶随机差分方程\[X_t=a_1X_{t-1}+a_2X_{t-2}+...+a_pX_{t-p}+\epsilon_t\]且满足如下条件

  (1){\(\epsilon_t\)}为白噪声序列

  (2)\(a_p\neq0\)且\(E{X_s}{\epsilon_t=0},s<t\)

称该模型为\(p\)阶自回归模型，记为\(AR(p)\)，非负整数\(p\)称为自回归阶数，实系数\(a_1,a_2,...a_p\)称为自回归序列系数。

3.2.2 MA模型

  对于自相关系数q阶截尾，偏自相关系数拖尾的数据，我们采用MA模型。设\(\left\{\epsilon_t\right\}\)为一白噪声序列，称以下模型为q阶滑动平均模型，简记为MA(q)模型（Moving Average),\[X_t=\epsilon_t-\beta_1\epsilon_{t-1}-...-\beta_q\epsilon_{t-q},t \in Z\]称序列\(\left\{\epsilon_t\right\}\)为q阶滑动平均序列，简称MA(q)序列。

注：白噪声序列通俗来说是零均值、常方差的稳定随机序列。

3.2.3 ARMA模型

  ARMA模型可以看成是由自回归模型（简称AR模型）与滑动平均模型（简称MA模型）为基础“混合”构成。

  对于偏自相关系数和自相关系数都拖尾的数据，我们采用ARMA模型。我们称以下模型为p阶自回归与q阶滑动平均混合模型,以后简记为ARMA(p,q)模型，即\[X_t=a_1X_{t-1}+a_2X_{t-2}+...+a_pX_{t-p}+\epsilon_t-\beta_1\epsilon_{t-1}-...-\beta_q\epsilon_{t-q}\]式中\(\left\{\epsilon_t\right\}\)为一白噪声序列，而且\(E(X_t\epsilon_t)=0\),对一切\(s<t\)成立。记\[a(u)=1-a_1u-...-a_pu^p\]\[b(u)=1-\beta_1u-...-\beta_qu^q\]\(a(u)\)与\(b(u)\)没有公共根。

4 非平稳序列

  前面介绍了对平稳时间序列的分析方法。实际上，在自然界中绝大部分序列都是非平稳的，因而对非平稳序列的分析更普遍、更重要，人们创造的分析方法也很多。这些分析方法总的来说可以分为两大类，确定性时序分析和随机时序分析，因为任何一个序列的波动都可以视为同时受到了确定性影响和随机性影响。平稳序列要求这两方的影响都是稳定的，而非平稳序列产生的机理就在于它所受到的这两方面的影响至少有一方面是不稳定的。

4.1 确定性时序分析

  由确定性因素导致的非平稳，通常显示出非常明显的规律性，这种规律性信息通常也比较容易提取，传统的时序分析方法一般会把分析的重点放在确定性信息的提取上。最常用的确定性分析方法是确定因素分解法。主要是说，序列的各种变化都可以归纳成四大类因素的综合影响。

  • 长期趋势（Trend）:时间序列随时间而逐渐增加或减少的变化趋势
  • 季节变动（Season）:时间序列在一年中的固定时间内，呈现出周期性重复变动
  • 循环变动（Circle）:沿着趋势线随钟摆般循环变动
  • 随机波动（Immediate）:受随机因素影响所引起的变动，没有规律

  而常见的时间序列主要包括3个部分：整体趋势、季节性和噪声。

  在确定性时序分析中，我们介绍一种常用的预测方法：平滑法。平滑法是趋势分析的一种方法，它是利用修匀技术，削弱短期随机波动对序列的影响，使序列平滑化，从而显示出变化的规律。它具有调节灵活、计算简便的特征，广泛应用于计量经济、人口研究等诸多领域。根据所用的平滑技术的不同，平滑法常用的有移动平均法和指数平滑法。

4.1.1 移动平均法

  移动平均法的基本思想是对于一个时间序列\(\left\{X_t\right\}\)，我们可以假定在一个比较短的时间间隔里，序列的取值是比较稳定的，它们之间的差异主要是由随机波动造成的。根据这种假定，我们可以用一定的时间间隔内的平均值作为某一期的估计值。

  • \(n\)期中心移动平均\[\tilde{x_t}=\begin{cases}\frac{1}{n}(x_{t-\frac{n-1}{2}}+x_{t-\frac{n-1}{2}+1}+...+x_t+...+x_{t+\frac{n-1}{2}-1}+x_{t+\frac{n-1}{2}}),n为奇数\\\frac{1}{n}(\frac{1}{2}x_{t-\frac{n}{2}}+x_{t-\frac{n}{2}+1}+...+x_t+...+x_{t+\frac{n}{2}-1}+\frac{1}{2}x_{t+\frac{n}{2}}),n为偶数\end{cases}\]

  股市中普遍使用的均线，帮助交易者确认股票或大盘的现有趋势，并预测其未来短期、中期或长期的走势，实际上，这些均线也就是中心移动平均线。

  • \(n\)期移动平均\[\tilde{x_t}=\frac{1}{n}(x_t+x_{t-1}+...+x_{t-n+1})\]

  移动平均的期数对原序列的修匀效果影响很大，如果想得到长期趋势，就应该做期数较大的移动平均；如果想密切管制序列的短期趋势，就应该做期数较小的移动平均。

4.1.2 指数平滑法

  移动平均法实际上就是用一个简单的加权平均数作为某一期趋势的估计值。以\(n\)期移动平均为例，\[\tilde{x_t}=\frac{1}{n}(x_t+x_{t-1}+...+x_{t-n+1})\]  相当于用近\(n\)期的加权平均数作为最后一期趋势的估计值，它们的权重都取为\(\frac{1}{n}\)，实际上也就是假定无论时间的远近，这\(n\)期的观察值对最后一期的影响都是一样的。但在实际生活中，我们会发现对大多数随机事件而言，一般都是近期的结果对现在影响会大些，远期的结果对现在的影响会小些。为了更好地反映这种影响作用，我们将考虑到时间间隔对事件发展的影响，各期权重随时间间隔的增大而呈指数衰减。这就是指数平滑法的基本思想。

  • 简单指数平滑

\[\tilde x_t=\alpha x_t+\alpha(1-\alpha)x_{t-1}+\alpha(1-\alpha)^2x_{t-2}+...\]

  式中，\(\alpha\)为平滑系数，它满足\(0<\alpha<1\)，平滑系数由研究人员根据经验给出，一般对于变化缓慢的序列，常取较小值，对于变化迅速的序列，常取较大的值。经验表明，平滑系数介于0.05~0.3之间的修匀效果比较好。从理论上我们可以证明，使用指数平滑法，预测任意l期的预测值都是常数，所以使用指数平滑法最好只做1期预测。

  • Holt两参数指数平滑

  该方法适用于对含有线性趋势的序列进行修匀。它的基本思想是假定序列有一个比较固定的线性趋势———每期都递增\(r\)或递减\(r\)，那么第\(t\)期的估计值就应该等于第\(t-1\)期的观察值加上每期固定的趋势变动值，即\(\hat x_t=x_{t-1}+r\)，但是由于随机因素的影响，使得每期的递增或递减值不会恒定为\(r\)，它会随时间变化上下波动，所以趋势序列实际上是一个随机序列{\(r_t\)} 。

  考虑用第\(t\)期的观察值和第\(t\)期的估计值的加权平均数作为第\(t\)期的修匀值：\[\tilde x_t=\alpha x_t+(1-\alpha)\hat x_t\]，为了让修匀序列{\(\tilde x_t\)}更平滑，我们对{\(r_t\)}也进行一次修匀处理：\[r_t=\gamma(\tilde x_t-\tilde x_{t-1})+(1-\gamma)r_{t-1}\]

  \(\alpha,\gamma\)是两个平滑参数，参数的选择原则也和简单指数平滑法一样，平滑序列的初始值最简单的是指定\(\tilde x_0=x_1\)，趋势序列初始值最简单的指定是用任意指定的区间的平均趋势作为趋势初始值：\(r_0=\frac{x_{n+1}-x_1}{n}\)。

4.2 随机时序分析

  非平稳序列的确定性因素分解方法原理简单、操作简便、易于解释，在宏观经济管理与预测领域有着广泛的应用。但是这种方法只能提取出强劲的确定性信息，对随机性信息浪费严重，导致模型拟合精度通常不够理想。随机时序分析方法的发展就是为了弥补确定性因素分解方法的不足。在这里我们介绍一种常用的随机时序建模方法：ARIMA模型。

4.2.1 ARIMA模型

  具有如下结构的模型称为ARIMA（求和自回归移动平均）模型，简记为ARIMA(p,d,q)模型：\[\begin{cases}\Phi(B)\nabla^dx_t=\Theta(B)\epsilon_t\\E(\epsilon_t)=0,Var(\epsilon_t)=\sigma_{\epsilon}^2,E(\epsilon_t\epsilon_s)=0,s\ne t\\Ex_s\epsilon_t=0,\forall s<t\end{cases}\]   式中，\(\nabla^d=(1-B)^d\);\(\Phi(B)=1-\phi_1B-...-\phi B^p\)为平稳可逆ARMA(p,q)模型的自回归系数多项式；\(\Theta(B)=1-\theta_1B-...-\theta_qB^q\)，为平稳可逆ARMA(p,q)模型的移动平滑系数多项式。

  可以看出，ARIMA模型的实质就是差分运算与ARMA模型的组合。这说明任何非平稳序列如果能通过适当阶数的差分实现差分后平稳，就可以对差分后序列进行ARMA模型拟合了，而ARMA模型的分析方法非常成熟，这意味着对差分平稳序列的分析也将是非常简单可靠的。

5 R的实现

5.1 相关软件包

TTR软件包,该软件包包含许多受欢迎的技术分析功能，例如我们今天要用到的简单移动平均函数就在这个包中。

forecast软件包,包含用于时间序列数据的预测函数

5.2 核心函数介绍

ts()函数

  ts()函数把数据转换成时间序列数据，调用方法为

ts(data = NA, start = 1, end = numeric(), frequency = 1,
   deltat = 1, ts.eps = getOption("ts.eps"), class = , names = )
  • data是时间序列观测值的一个数值向量或矩阵，数据框将通过data.matrix被强制转换成数值向量；
  • start指定收集数据的第一年及第一个间隔期，例如我们想指定第一个时间为2017年二季度，则设置start=c(2017,2)；
  • end指定时间序列的终止时间点；
  • 有些时间序列数据是间隔小于一年的数据，如月度或季度数据，参数frequency就用来指定数据在一年中的频数，例如月度数据设置为frequency=12;
  • delta表示连续观测值周期的分数，与frequency的值互为倒数。

is.ts()函数

  判断数据是否是时间序列数据，返回值是逻辑值，若值为TRUE，则表示是时间序列数据

plot.ts()函数

  绘制时间序列图

decompose()和stl()函数

  估计时间序列的趋势、季节性和不规则部分，函数调用方法分别为

decompose(x, type = c("additive", "multiplicative"), filter = NULL)
  • x是时间序列对象；
  • type指定季节分解的形式，是相加还是相乘；
  • filter是滤波系数，用T表示趋势，S表示周期，e表示噪声。
stl(x,s.window,s.degree=0,...)
  • x是时间序列对象；
  • s.window为字符串“periodic”或“”Loess“方法提取季节的跨度，应为奇数值；
  • s.degree是局部多项式拟合季节性提取的程度，取值为1或0。

diff()函数

  对数据进行差分处理，调用方法是

diff(x,lag=1,differences=1,...)
  • x是待差分的数值向量或矩阵；
  • lag是指定差分的阶数

SMA()函数

  是实现简单移动平均的函数，调用方法是

SMA(x,n=10,...)
  • n是指定移动平均的跨度，跨度越大，估计得到的趋势部分越平滑。

auto.arima()函数

  是R语言提供的自动拟合ARIMA模型、AR模型、MA模型、ARMA模型的函数，拟合哪种模型根据某种准则确定的，准则可以是AIC准则、BIC准则，函数的调用方法是

auto.arima(x,ic=c("aicc", "aic", "bic"))

5.3 相关数据集

  我们使用R的内置数据集：Box&Jenkins航空公司1949-1960年每月国际航线乘客数，首先来看一下数据集。

data("AirPassengers")#获取数据集
data=AirPassengers#为简化变量名，将数据集命名为data
is.ts(data)#判断是否为时间序列数据
## [1] TRUE
plot.ts(data)#绘制时间序列图

6 应用实例

  探索适合该数据集的时间序列分析模型。首先可以看出数据具有明显的周期性和趋势性，用decompose()函数来估计时间序列的趋势、季节性和不规则部分。

co=decompose(data)
options(dights=3)#显示小数点后3位有效数字
plot(co)#画出估计的时间序列趋势、季节性和不规则部分

  图中包含了4部分，从上到下分别是：原始的时间序列观测值；估计出的趋势部分；估计出的季节性部分；估计出的随机波动。从这几张图中可以看得非常明显，整体呈现增长的趋势，其中，1954-1958年间增长迅速；周期性十分规则。其实除了decompose()函数可以用于分解周期性时间序列外，还有stl()函数调用方法前面有介绍过，呈现结果大同小异，大家可以动手试试。

  用简单移动平均对序列修匀，减小随机波动的影响，在R中实现简单移动平均的函数是SMA()函数，该函数在TTR程序包中。

#install.packages("TTR"")
library(TTR)#加载软件包
SMA.data=SMA(data,n=3)#调用SMA函数
plot(SMA.data,n=3)#绘制时间序列图像

  使用跨度为3的简单易懂平均法后，可以看到图像确实平滑了许多。我们可以看到乘客数有一个比较固定的增长趋势，因此，可以用Holt两参数平滑法进行预测短期预测

#install.packages(forecast)#下载forecast软件包
library(forecast)#加载软件包
data.pre=HoltWinters(data,gamma=FALSE)
data.pre#显示预测的各项结果
## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = data, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 1
##  beta : 0.003218516
##  gamma: FALSE
## 
## Coefficients:
##         [,1]
## a 432.000000
## b   4.597605
plot(data.pre)#绘图对比Holt平滑对原始数据的预测效果

  可以看出，原始数据和预测值趋势相似，差别不大，预测值似乎会早一步比原始数据出现。接下来用forecast包中的forecast()函数来预测10年后的情况。

data.pre2=forecast(data,h=24)#预测两年，所以是h取24个月份
plot(data.pre2)

  阴影部分显示预测结果，中间较粗的曲线是24个时点的乘客数的点估计值，深色阴影是置信水平为80%的预测区间，浅色阴影是90%的预测区间。有兴趣的同学可以查查1961-1962年间的真实数据，看看我们的预测结果跟真实情况有多大出入呢。

  接下来我们用随机时序分析方法对数据进一步探索。ARIMA模型是将非平稳序列化为平稳序列的操作，具体是通过差分的方法，ARIMA(p,d,q)中的d就是差分阶数，我们来看看差分的效果，先进行一阶差分。

data.diff=diff(data)#做差分
plot.ts(data.diff)#绘制时间序列图

  从图形中可以看出来，序列波动的范围是越来越大的，说明序列的方差是越来越大的，我们通常会用取对数的方式对序列进行去异方差性处理。

data1=log(data)#取对数
data1.diff=diff(data1)#一阶差分
plot(data1.diff)#绘制差分后时间序列图

  通过自相关系数图观察序列是否平稳

acf(data1.diff,lag.max = 10)#绘制自相关系数图，lag.max是最大滞后阶数

  虽然时间序列图看似随机波动状态，没有明显趋势周期，但是自相关系数图表明，10阶之内还有一些阶数并没有落在置信区间内（蓝线表示置信区间），因此我们尝试再一次差分。

data2=diff(data.diff)
acf(data2,lag.max = 10)#绘制自相关系数图

  10阶以内的自相关系数基本都落在置信区间内，可以近似认为序列平稳.接下来，就可以用平稳序列建模的方法进行分析了.对ARIMA(p,d,q)中阶数p和q的判断，R提供了一个快捷方式，通过auto.arima()来发现恰当的ARIMA模型

auto.arima(data2,ic="aic")#使用aic准则
## Series: data2 
## ARIMA(5,0,0)(0,1,0)[12] 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ar5
##       -1.0906  -0.8005  -0.6538  -0.5638  -0.2573
## s.e.   0.0845   0.1184   0.1240   0.1178   0.0858
## 
## sigma^2 estimated as 167.8:  log likelihood=-515.72
## AIC=1043.44   AICc=1044.12   BIC=1060.64

  最终，R为我们选择了ARIMA(0,1,1)模型，说明R认为我们应该在之前的基础上再一次差分。所以我们的结论是，对数据进行对数变换后，再进行三次差分，对差分后的序列拟合MA（移动平均）模型。

7 本章汇总

名称 类别 功能
TTR 软件包 提供函数SMA()
forecast 软件包 时间序列中的通用预测函数
AirPassengers 数据集 某航空公司每月国际航线乘客数据
is.ts() 函数 判断数据是否为时间序列数据
plot.ts() 函数 绘制时间序列图
decompose() 函数 用于分解时间序列
options() 函数 控制有效数字输出
stl() 函数 与decompose()函数作用相同
SMA() 函数 实现简单移动平均
HoltWinters() 函数 用于指数平滑法修匀和预测
diff 函数 差分操作
log() 函数 取对数操作
acf() 函数 画自相关系数图
auto.arima() 函数 自动拟合ARIMA模型的函数

8 参考文献

[1]安鸿志.时间序列分析[M],华东师范大学出版社,1992.

[2]王燕.应用时间序列分析[M],中国人民大学出版社,2012.

[3]李诗羽.数据分析：R语言实战[M],电子工业出版社,2014

源文件下载 在线调试
返回： 常用统计学, 时间序列
Prev Post
Next Post

2025@北京大学古生物形态学定量分析实践平台 All Rights Reserved
北京大学地球与空间科学学院
北京大学地球科学国家级实验教学示范中心
北京大学地球科学国家级虚拟仿真实验教学中心