【译】时间序列建模完整教程（R语言）

对于企业时间是最重要的因素，然而绝大多数公司很难跟上时间的脚步。但是随着技术的发展，出现了很多有效的方法，能够让我们预测未来。不要担心，本文并不会讨论时间机器，讨论的都是很实用的东西。

本文将要讨论关于预测的方法。有一种预测是跟时间相关的，而这种处理与时间相关数据的方法叫做时间序列模型。其中一种处理基于时间的数据的方法是时间序列建模。 顾名思义，它涉及到按时(年、日、小时、分钟)的数据，以获得隐藏的见解，从而做出明智的决策。

当我们处理时序序列数据的时候，时间序列模型是非常有用的模型。大多数公司都是基于时间序列数据来分析第二年的销售量，网站流量，竞争地位和更多的东西。然而很多人并不了解的时间序列分析这个领域。 所以，如果你不了解时间序列模型。这篇文章将会想你介绍时间序列模型的处理步骤以及它的相关技术。

本文包含的内容如下所示：

• 时间序列模型介绍
• 使用R语言来探索时间序列数据
• 介绍ARMA时间序列模型
• ARIMA时间序列模型的框架与应用

什么是时间序列建模

让我们从基础开始。这包括平稳序列、随机行走、 Rho 系数、 Stationarity 的 Dickey Fuller 平稳性检验。如果这些术语已经吓到你了，不要担心——它们很快就会变得清晰起来，我敢打赌你会在我解释的时候开始享受这个主题。

平稳序列

判断一个序列是不是平稳序列有三个评判标准：

1、均值，是与时间t 无关的常数。下图（左）满足平稳序列的条件，下图（右）很明显具有时间依赖。

2、方差，是与时间t 无关的常数。这个特性叫做方差齐性。下图显示了什么是方差对齐，什么不是方差对齐。（注意右手边图中的不同分布）

3、协方差，只与时期间隔k有关，与时间t 无关的常数。如下图（右），可以注意到随着时间的增加，曲线变得越来越近。因此红色序列的协方差并不是恒定的。

我们为什么要关心时间序列的“平稳性”？

除非你的时间序列是平稳的，否则不能建立一个时间序列模型。在很多案例中时间平稳条件常常是不满足的，所以首先要做的就是让时间序列变得平稳，然后尝试使用随机模型预测这个时间序列。有很多方法来平稳数据，比如消除长期趋势，差分化。

随机漫步

这是时间序列中最基本的概念。 你可能很了解这个概念。 但是，我发现在这个行业里，很多人把随机漫步理解为平稳过程。 在这一节中，我会使用一些数学工具，帮助理解这个概念。我们先看一个例子。

例如: 想象一个女孩在一个巨大的棋盘上随意移动。 在这种情况下，女孩的下一个位置只取决于最后一个位置。

 

这里的 代表这这个时间点随机干扰项。这个就是女孩在每一个时间点带来的随机性。

现在我们递归所有x时间点，最后我们将得到下面的等式：

 

现在，让我们尝试验证一下随机游走的平稳性假设：

1、是否均值为常数？

 

我们知道由于随机过程的随机干扰项的期望值为0.到目前为止：E[X(t)] = E[X(0)] = 常数

2、方差是恒定的吗？

 

 

因此，我们推断，随机漫步不是一个平稳过程，因为它具有时变方差。另外，如果我们检查协方差，我们看到协方差也依赖于时间。

Rho系数

我们已经知道一个随机游走是一个非平稳的过程。让我们在方程中引入一个新的系数，看看我们是否能制定一个检查平稳性的公式。

引入系数: Rho

 

现在，我们将改变Rho看看我们可不可以让这个序列变的平稳。这里我们将从视觉上进行确定，而不是做任何测试来检查平稳性。

让我们从一个Rho=0的完全平稳序列开始。这里是时间序列的图：

将Rho的值增加到0.5，我们将会得到如下图：

你可能会注意到，周期变长了，但基本上似乎没有一个严重的违反平稳性假设。现在让我们采取更极端的情况下ρ= 0.9

我们仍然看到，在一定的时间间隔后，从极端值返回到零。这一系列也不违反非平稳性明显。现在，让我们用ρ= 1随机游走看看

这显然是违反平稳条件。是什么使rho=1变得这么特殊的呢？，这种情况并不满足平稳性测试？我们来找找这个数学的原因

公式 的期望为：

这个公式很有意义。下一个X(或者时间点t)被拉到Rho*上一个x的值。

例如，如果x(t–1)=1，E[X(T)] = 0.5（Rho=0.5）。现在，如果从零移动到任何方向下一步想要期望为0。唯一可以让期望变得更大的就是错误率。当Rho变成1呢？下一步没有任何可能下降。

Dickey Fuller平稳性检验

这里学习的最后一个知识点是Dickey Fuller检验。。在统计学里，Dickey-Fuller检验是测试一个自回归模型是否存在单位根。这里根据上面Rho系数有一个调整，将公式转换为Dickey-Fuller检验

 

我们要测试如果Rho–1=0是否差异显著。如果零假设不成立，我们将得到一个平稳时间序列。 平稳性测试和将一个序列转换为平稳性序列是时间序列模型中最重要的部分。因此需要记住本节提到的所有概念方便进入下一节。

使用R语言探索时间序列

本节我们将学习如何使用R处理时间序列。这里我们只是探索时间序列，并不会建立时间序列模型。 本节使用的数据是R中的内置数据：AirPassengers。这个数据集是1949-1960年每个月国际航空的乘客数量的数据。

加载数据集

下面的代码将帮助我们载入数据集并且能够看到一些数据集的宏观指标。

> data(AirPassengers) #载入数据
> class(AirPassengers) #查看AirPassengers数据类型，这里的"ts"代表的是时间序列数据
[1] "ts"
> start(AirPassengers) #查看开始时间
[1] 1949    1
> end(AirPassengers) #查看结束时间
[1] 1960   12
> frequency(AirPassengers) #查看时间概率的频率，这里是一年12个月
[1] 12
> summary(AirPassengers) #查看数据的一些宏观指标
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  104.0   180.0   265.5   280.3   360.5   622.0

查看详细数据

> plot(AirPassengers) #绘制出时间序列
> abline(reg=lm(AirPassengers~time(AirPassengers))) #拟合一条直线

> cycle(AirPassengers)
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949   1   2   3   4   5   6   7   8   9  10  11  12
1950   1   2   3   4   5   6   7   8   9  10  11  12
1951   1   2   3   4   5   6   7   8   9  10  11  12
1952   1   2   3   4   5   6   7   8   9  10  11  12
1953   1   2   3   4   5   6   7   8   9  10  11  12
1954   1   2   3   4   5   6   7   8   9  10  11  12
1955   1   2   3   4   5   6   7   8   9  10  11  12
1956   1   2   3   4   5   6   7   8   9  10  11  12
1957   1   2   3   4   5   6   7   8   9  10  11  12
1958   1   2   3   4   5   6   7   8   9  10  11  12
1959   1   2   3   4   5   6   7   8   9  10  11  12
1960   1   2   3   4   5   6   7   8   9  10  11  12
> plot(aggregate(AirPassengers,FUN=mean)) #绘制周期内（每年）的平均值
> boxplot(AirPassengers~cycle(AirPassengers)) #绘制箱形图，了解周期内的数据变化（最大值、最小值、中位数、及上下四分位数。）

重要推论

• 趋势显示旅客的数量每年都在增加
• 七八月的均值和方差比其他月份要高很多
• 每个月的平均值并不相同，但是方差差异很小。因此，可以看出具有很强的周期性。一个周期为12个月或更少。

在时间序列模型中，探索数据是最重要的一步——如果没有这一步，你将不知道这个序列是不是平稳序列。就像这个例子一样，我们已经知道了很多关于这个模型的很多细节。

ARMA 时间序列建模

ARMA也叫自回归移动平均混合模型。ARMA模型经常在时间序列中使用。在ARMA模型中，AR代表自回归，MA代表移动平均。如果这些术语听起来很复杂，不用担心-下面将会用几分钟的时间简单介绍这些概念。

在开始之前，你首先要记住，AR或者MA并不是应用在非平稳序列上的。 在实际应用中可能会得到一个非平稳序列，你首先要做的就是将这个序列变成平稳序列（通过差分化/转换），然后选择可以使用的时间序列模型。

首先，本文将介绍分开介绍两个模型（AR&MA）。接下来我们看一看这些模型的特点。

自回归时间序列模型（Auto-Regressive Time Series Model）

让我们用下面的例子来理解 AR 模型:一个国家当前的GDP称x(t)取决于去年的GDP，即x(t-1)。一个国家在一个财政年度的GDP取决于前一年建立的制造厂/服务，以及本年度新成立工厂/服务。但是GDP的主要组成部分是前一部分，因此，我们可以把 GDP 的等式写成:

 

这个方程称为AR公式。公式表示下一个点完全依赖与前面一个点。α是我们寻求的系数，以便最小化误差函数。注意，x(t-1)同样依赖 x(t-2)。 因此，波动在未来会逐渐消退。

例如，假设x(t)是一个城市在某一天出售的果汁瓶的数量。 在冬天，很少有供应商购买果汁瓶。突然，在一个特定的日子里，温度上升，果汁的需求猛增到1000个。然而，几天后，气候又变得寒冷起来。但是，知道人们在炎热的日子里已经习惯了喝果汁，在寒冷的日子里仍有50% 的人还在喝果汁。接下来的几天，这个比例下降到25% (50%的50%)，然后在数天之后逐渐变小。下面的图表说明了AR 系列的惯性特性：

移动平均时间序列模型（Moving Average Time Series Model）

让我们以另一个例子来理解移动平均时间序列模型：一个公司生产了一批类型在市场上很容易买到的包，在竞争激烈的市场。这个包的销量在很多天都是0，于是，在某一天，这家公司做了一个实验，生产了一个不同类型的包，这种包在市场上任何地方找不到，因此，他可以一下子卖掉了1000个。这个包的需求特别高，很快库存快要完了。这天结束了还有100个包没卖掉。我们把这个误差成为时间点误差。接下来的几天仍有几个客户购买这种包。随着时间的推移，这个袋子已经失去了它的吸引力。下面通过一个简单的公式来描述这个场景：

 

如果我们尝试绘制这个图表，它会看起来像这样：

注意到 MA 和 AR 模型的区别了吗？ 在 MA 模型中，噪声/冲击随着时间的推移而迅速消失。 AR模型对冲击具有持久的影响。

AR模型与MA模型的不同

AR和MA 模型之间的主要区别是基于不同时间点的时间序列对象之间的相关性。x(t)与x(t-n)的相关性总为0，这直接源于 x(t)和 x(t-n)之间的协方差对于MA模型来说是零。然而，AM模型中x(t)与x(t-1)的相关性随着时间的推移变得越来越小。不管是否有 AR 模型或 MA 模型，这种差异都可以被利用。

利用ACF和PACF绘图

一旦我们得到平稳时间序列，我们必须回答两个主要问题：

• Q1：是AR 还是MA 过程？
• Q2：我们需要使用哪种顺序的AR 或MA 过程？

解决这些问题的诀窍可在上一节中找到。你没注意到吗？

第一个问题可以用总相关图（也称为自相关函数/ACF）来回答。ACF是不同滞后函数之间的总关联图。例如，在GDP问题上，时间点T的GDP是x(t)。我们对x(t)与x(t-1)、x(t-2)等的相关性很感兴趣。现在让我们来反思一下上面所学的内容。在移动平均序列的滞后n中，x(t)与x(t-n-1)之间不存在任何相关。因此，总相关图在第n 次滞后处截止。因此，很容易找到一个MA系列的滞后。对于AR 序列，这种相关性将逐渐下降而没有任何截止值。那么，如果我们是AR 系列，我们该怎么办？

这是第二个窍门。如果发现每个滞后的部分相关，则在AR 级数的程度之后会被切断。例如，如果我们有AR系列，如果我们排除第一滞后x(t-1)的影响，则我们的第二滞后x(t-2)与x(t)无关。因此，偏相关函数（PACF）在第一滞后急剧下降。下面是一些例子来澄清你对这个概念的任何怀疑：

上面的蓝线显示出明显不同于零的值。显然，上面的图在第二滞后之后在PACF 曲线上有一个截止，这意味着这主要是AR过程。

显然，上面的图在第二滞后之后在ACF 曲线上有一个截止，这意味着这主要是MA过程。

到目前为止，我们已经讨论了如何识别类型的固定系列使用ACF 和PACF 图。现在，我将向大家介绍一个构建时间序列模型的全面框架。此外，我们还讨论了时间序列建模的实际应用。

ARIMA时间序列模型的框架与应用

到此，本文快速介绍了时间序列模型的基础概念、使用R探索时间序列和ARMA模型。现在我们将这些零散的东西组织起来，做一件很有趣的事情。

框架概述

下图的框架展示了如何一步一步的“做一个时间序列分析”

前三步我们在前文已经讨论了。尽管如此，这里还是需要简单说明一下：

第一步：时间序列可视化

在构建任何类型的时间序列模型之前，分析其趋势是至关重要的。我们感兴趣的细节与系列中的任何趋势、季节性或随机行为有关。我们已经在本系列的第二部分中涵盖了这一部分。

第二步：序列平稳

一旦我们知道了模式、趋势、周期。我们就可以检查序列是否平稳。Dicky-Fuller是一种很流行的检验方式。在第一部分意见介绍了这种检验方式。在这里还没有结束！如果发现序列是非平稳序列怎么办？

这里有三种比较常用的技术来让一个时间序列平稳。

1、消除趋势：这里我们简单的删除时间序列中的趋势成分。例如，我的时间序列的方程是：

 

我们只需删除括号中的部分，并为其余部分建立模型。

2、差分：这是常用的去除非平稳性的技术。这里我们是对序列的差分的结果建立模型而不是真正的序列。例如：

 

这种差异称为 AR(I)MA 中的积分部分。 现在，我们有三个参数：

p:AR
d:I
q:MA

3、季节性：季节性直接被纳入ARIMA模型中。下面的应用部分我们再讨论这个。

第三步：找到最优参数

使用ACF 和PACF 图可以找到参数p、d、q。这种方法的一个补充是，如果ACF 和PACF 逐渐减小，则表明我们需要使时间序列平稳并将值引入到“d”。

第四步：建立ARIMA模型

找到了这些参数，我们现在就可以尝试建立ARIMA模型了。从上一步找到的值可能只是一个近似估计的值，我们需要探索更多(p,d,q）的组合。最小的BIC和AIC的模型参数才是我们要的。我们也可以尝试一些具有季节性成分的模型。在这里，在ACF/PACF图中我们会注意到一些季节性的东西。

第五步：预测

一旦我们有了最后的ARIMA 模型，我们现在就可以预测未来的时间点。如果模型工作正常，我们也可以用可视化趋势来验证。

时间序列模型的应用

这里我们用前面的例子。使用这个时间序列做预测。我们建议你在进行下一步之前，先观察这个数据。

下图是这些年的乘客数的图。在往下看之前，观察这个图。

这里是我的观察：

• 乘客有着逐年增加的趋势。
• 这看起来有季节性，每一个周期不超过12个月。
• 数据的方差逐年增加。

在我们进行平稳性测试之前我们需要解决两个问题。第一，我们需要消除方差不齐。这里我们对这个序列做求对数。第二我们需要解决序列的趋势性。我们通过对时序序列做差分。现在，我们来检验最终序列的平稳性。

> adf.test(diff(log(AirPassengers)), alternative="stationary", k=0)
Error in adf.test(diff(log(AirPassengers)), alternative = "stationary",  : 
  could not find function "adf.test"

验证时发现没有安装”tseries” 包，解决方案：

> install.packages("tseries")
> library(tseries)

测试结果：

	Augmented Dickey-Fuller Test
 
data:  diff(log(AirPassengers))
Dickey-Fuller = -9.6003, Lag order = 0, p-value = 0.01
alternative hypothesis: stationary
 
Warning message:
In adf.test(diff(log(AirPassengers)), alternative = "stationary",  :
  p-value smaller than printed p-value

下一步是找到在ARIMA 模型中使用的正确参数。我们已经知道，d分量是1，因为我们需要1的差值使级数平稳。我们使用相关图来实现这一点。以下是ACF 系列图：

acf(log(AirPassengers))

在上面的图表中你看到了什么？很显然ACF下降的十分的慢，这就意味着乘客的数量并不是平稳的。我们在前面已经讨论了，我们现状准备在序列去对数后的差分上做回归，而不是直接在序列去对数后的数据熵差分。让我们看一下差分后的ACF和PACF曲线吧。

> acf(diff(log(AirPassengers)))

pacf(diff(log(AirPassengers)))

显然，ACF 曲线在第一次滞后之后就被切断了。因此，我们理解P 的值应该是0，因为ACF是被切断的曲线。Q 值为1 或2。经过几次迭代，我们发现（0，1，1）AS（p，d，q）是最小AIC 和BIC 的组合。

让我们建立一个 ARIMA 模型，并预测未来的10年。 此外，我们将尝试在 ARIMA 的配方中加入季节性成分。 然后，我们将把预测和训练数据一起显示出来。

> (fit <- arima(log(AirPassengers), c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12)))
 
Call:
arima(x = log(AirPassengers), order = c(0, 1, 1), seasonal = list(order = c(0, 
    1, 1), period = 12))
 
Coefficients:
          ma1     sma1
      -0.4018  -0.5569
s.e.   0.0896   0.0731
 
sigma^2 estimated as 0.001348:  log likelihood = 244.7,  aic = -483.4
> pred <- predict(fit, n.ahead = 10*12)
> ts.plot(AirPassengers,2.718^pred$pred, log = "y", lty = c(1,3))

原文链接： https://www.analyticsvidhya.com/blog/2015/12/complete-tutorial-time-series-modeling/