Skip to content

Latest commit

 

History

History
767 lines (485 loc) · 48 KB

2.md

File metadata and controls

767 lines (485 loc) · 48 KB
Raw

二、了解时间序列数据

在上一章中，我们介绍了时间序列分析的一般方法，该方法包括两个主要步骤：

  • 数据可视化以检查趋势，季节性和周期性模式的存在
  • 调整趋势和季节性以生成平稳序列

生成固定数据对于增强时间序列预测模型非常重要。 趋势，季节和周期性成分的推论将使我们面临不规则的波动，而不能仅以时间指数作为解释变量来建模。 因此，为了进一步改善预测，不规则波动被假定为独立于且分布均匀的（iid）观测值，并通过对时间指数以外的变量进行线性回归建模。

例如，房价可能同时呈现趋势和季节性（例如，季度）变化。 但是，在调整趋势和季节性之后剩余的残差实际上可能取决于外生变量，例如总建筑面积，建筑物中的楼层数等，这取决于采样数据中的特定实例。 因此，趋势和季节性调整以及外生变量模型将是对时间序列未来实例的更好预测。

将原始时间序列更改为 iid 观测值，或者换句话说，使时间序列平稳，是开发基于外生变量的线性回归模型的重要步骤。 这是因为存在完善的统计方法，例如中心极限定理，最小二乘法等，这些方法对于 iid 观测非常有效。

下面的流程图总结了前面各段中描述的时间序列分析方法。 在本章中，我们将在以下主题中介绍此方法的步骤 1、2 和 3：

  • 时间序列数据的高级处理和可视化
  • 统计假设检验以验证时间序列的平稳性
  • 时间序列分解以调整趋势和季节性

读者会发现本章讨论的数学概念是开发时间序列预测模型的基本构建块。

图 2.1：时间序列分析的通用方法

时间序列数据的高级处理和可视化

在许多情况下，原始时间序列需要转换为汇总统计信息。 例如，原始时间序列中的观测值可能每秒记录一次。 但是，为了执行任何有意义的分析，必须每分钟汇总一次数据。 这将需要在比原始数据中的细化时间索引更长的时间段内对观察值进行重采样。 对于每个较长的时间段，都会计算汇总统计信息，例如均值，中位数和方差。

时间序列数据预处理的另一个示例是计算数据中相似段上的聚合。 考虑由 X 公司制造的汽车的月度销售情况，该数据显示每月的季节性变化，因此给定年份的一个月内的销售情况类似于上一年度和下一年的同一月的销售情况。 为了突出这种季节性，我们必须从数据中删除长期趋势。 但是，让我们假设没有长期趋势或已经对其进行了调整。 我们有兴趣估算季节（在本例中为月）的统计数据，以确定季节之间（月之间）的变化。

例如，我们要求在 1 月，2 月等期间的平均销售额，以了解特定年份的平均销售额变化情况。 通过首先将原始时间序列分为 12 个细分（每个细分由一个月定义），然后汇总每个细分的数据，可以获得按季节划分的趋势（例如按月销售汽车）。

请注意，重采样和分组操作都将原始时间序列划分为不重叠的块，以找到聚合。 两种技术都将减少噪声并随后产生原始时间序列的平滑。

但是，有时需要时间序列的连续或连续汇总才能进行分析。 在连续时间段的窗口上计算聚合的技术给出了移动或滚动的聚合。 例如，汽车销售数据的季度移动平均值将是四个月窗口内的平均值，该窗口每次均会移动一个月。 通过移动或滚动计算窗口，可以生成移动平均值。

以下三个子节中通过示例演示了这些技术：

  • 重采样时间序列数据
  • 执行group-by
  • 计算移动统计

我们将使用 pandas 数据处理 API 来实现这些技术。

重采样时间序列数据

重新映射的技术使用了从 1975 年 1 月 1 日 st 到 1975 年 1 月 17 日之间每两小时获取的化学浓度读数的时间序列来说明。该数据集已从这个页面下载，也可以在本书的 GitHub 存储库的 datasets 文件夹中找到。

我们首先导入运行此示例所需的软件包：

from __future__ import print_function 
import os 
import pandas as pd 
import numpy as np 
%matplotlib inline 
from matplotlib import pyplot as plt

然后，我们将工作目录设置如下：

os.chdir('D:/Practical Time Series')

接下来，从pandas.DataFrame中的 CSV 文件中读取数据，并显示形状和DataFrame的前 10 行：

df = pd.read_csv('datasets/chemical-concentration-readings.csv') 
print('Shape of the dataset:', df.shape) 
df.head(10)

上面的代码返回以下输出：

Shape of the dataset: (197, 2)

| | 时间戳记 | 化学浓。 | | 0 | 1975-01-01 00:00:00 | 17.0 | | 1 | 1975-01-01 02:00:00 | 16.6 | | 2 | 1975-01-01 04:00:00 | 16.3 | | 3 | 1975-01-01 06:00:00 | 16.1 | | 4 | 1975-01-01 08:00:00 | 17.1 | | 5 | 1975-01-01 10:00:00 | 16.9 | | 6 | 1975-01-01 12:00:00 | 16.8 | | 7 | 1975-01-01 14:00:00 | 17.4 | | 8 | 1975-01-01 16:00:00 | 17.1 | | 9 | 1975-01-01 18:00:00 | 17.0 |

我们将通过在第二列上应用resample和mean函数将原始时间序列的每两小时观测值转换为每日平均值。 resample函数要求DataFrame的行索引为numpy.datetime64类型的时间戳。 因此，我们将行索引从整数（如上表所示）更改为datetime_rowid，这是numpy.datetime64对象的pandas.Series。 通过使用pd.to datetime实用程序功能从Timestamp列中生成numpy.datetime64对象。 以下代码显示了如何按行重新索引：

datetime_rowid = df['Timestamp'].map(lambda t: pd.to_datetime(t, format='%Y-%m-%d %H:%M:%S')) 
df.index = datetime_rowid 
df.head(10)

最后一行返回重新索引的DataFrame的前十行，如下所示。 请注意，新行索引的类型为numpy.datetime64。

| 时间戳记 | 时间戳记 | 化学浓。 | | | | | | 1975-01-01 00:00:00 | 1975-01-01 00:00:00 | 17.0 | | 1975-01-01 02:00:00 | 1975-01-01 02:00:00 | 16.6 | | 1975-01-01 04:00:00 | 1975-01-01 04:00:00 | 16.3 | | 1975-01-01 06:00:00 | 1975-01-01 06:00:00 | 16.1 | | 1975-01-01 08:00:00 | 1975-01-01 08:00:00 | 17.1 | | 1975-01-01 10:00:00 | 1975-01-01 10:00:00 | 16.9 | | 1975-01-01 12:00:00 | 1975-01-01 12:00:00 | 16.8 | | 1975-01-01 14:00:00 | 1975-01-01 14:00:00 | 17.4 | | 1975-01-01 16:00:00 | 1975-01-01 16:00:00 | 17.1 | | 1975-01-01 18:00:00 | 1975-01-01 18:00:00 | 17.0 |

现在我们准备在Chemical conc.列上应用resample和mean函数：

daily = df['Chemical conc.'].resample('D') 
daily_mean = daily.mean()

请注意，我们已将参数D传递给resample函数以生成每日平均值。 对于每月和每年的汇总，我们需要将M和Y传递给resample函数。

最后，原始平均值和每日平均值绘制在下图中，该图显示了后者的平滑效果：

图 2.2：每两小时读数的时间序列和化学浓度的每日平均值

生成上述图形并将其另存为 PNG 文件的代码如下：

fig = plt.figure(figsize=(5.5, 5.5)) 
ax = fig.add_subplot(1,1,1) 

df['Chemical conc.'].plot(ax=ax, color='b') 
daily_mean.plot(ax=ax, color='r') 

ax.set_title('Bi-hourly reading (blue) & Daily Mean (red)') 
ax.set_xlabel('Day in Jan 1975') 
ax.set_ylabel('Chemical concentration') 
plt.savefig('plots/ch2/B07887_02_17.png', format='png', dpi=300)

分组聚合

为了演示分组聚合，我们将使用美国得克萨斯州费舍尔河平均日温度的时间序列。 该时间序列具有从 1988 年 1 月 1 日 st 到 1991 年 12 月 31 日 st 之间的观测值。该数据集已从这个页面下载。 也可以在本书的 GitHub 存储库的 datasets 文件夹中找到。

我们首先阅读并绘制原始时间序列，如下所示：

图 2.3：美国德克萨斯州费舍尔河每日平均温度的时间序列

原始时间序列似乎具有每年重复的月度模式，可以通过计算月度平均值进行验证。 这是通过将数据分组为 12 个月，然后计算每个月的平均值来完成的。 以下代码段显示了执行此操作的代码。 我们首先在DataFrame中添加Month_Year列：

df['Month_Year'] = df.index.map(lambda d: d.strftime('%m-%Y'))

接下来，将Mean temparature列相对于新添加的Month_Year列进行分组，并计算按月的均值，中位数和标准差：

monthly_stats = df.groupby(by='Month_Year')['Mean temparature'].aggregate([np.mean, np.median,                            np.std]) 
monthly_stats.reset_index(inplace=True) 
monthly_stats.head(10)

下表显示按月汇总：

| | 月 _ 年 | 均值 | 中位数 | std | | 0 | 01-1988 | -22.137097 | -23.00 | 5.260640 | | 1 | 01-1989 | -17.129032 | -18.00 | 8.250725 | | 2 | 01-1990 | -15.112903 | -12.00 | 6.606764 | | 3 | 01-1991 | -23.038710 | -24.50 | 7.095570 | | 4 | 02-1988 | -19.025862 | -19.50 | 8.598522 | | 5 | 02-1989 | -19.267857 | -19.25 | 8.092042 | | 6 | 02-1990 | -17.482143 | -16.50 | 8.018477 | | 7 | 02-1991 | -10.967857 | -12.15 | 8.220753 | | 8 | 03-1988 | -8.258065 | -9.25 | 5.341459 | | 9 | 03-1989 | -12.508065 | -9.50 | 8.289925 |

请注意，上表中的行不是按Month_Year的升序排列。 因此，需要对其进行记录。 这是通过创建两个新列-Month和Year，然后以Year的升序排序，然后以Month的升序排序来完成的：

monthly_stats['Year'] = monthly_stats['Month_Year']\ 
                                      .map(lambda m: pd.to_datetime(m, format='%m-%Y').strftime('%Y')) 
monthly_stats['Month'] = monthly_stats['Month_Year']\ 
                        .map(lambda m: pd.to_datetime(m, format='%m-%Y').strftime('%m')) 
monthly_stats.sort_values(by=['Year', 'Month'], inplace=True) 
monthly_stats.head(10)

此处给出了已排序的DataFrame的前 10 行：

| | 月 _ 年 | 均值 | 中位数 | std | 年 | 月 | | 0 | 01-1988 | -22.137097 | -23.000 | 5.260640 | 1988 | 01 | | 4 | 02-1988 | -19.025862 | -19.500 | 8.598522 | 1988 | 02 | | 8 | 03-1988 | -8.258065 | -9.250 | 5.341459 | 1988 | 03 | | 12 | 04-1988 | 2.641667 | 1.875 | 5.057720 | 1988 | 04 | | 16 | 05-1988 | 11.290323 | 11.000 | 6.254364 | 1988 | 05 | | 20 | 06-1988 | 19.291667 | 19.000 | 3.909032 | 1988 | 06 | | 24 | 07-1988 | 19.048387 | 18.500 | 3.073692 | 1988 | 07 | | 28 | 08-1988 | 17.379032 | 18.000 | 3.183205 | 1988 | 08 | | 32 | 09-1988 | 10.675000 | 10.750 | 3.880294 | 1988 | 09 | | 36 | 10-1988 | 2.467742 | 3.000 | 6.697245 | 1988 | 10 |

下图绘制了月度总计，突出显示了原始数据中存在的按月季节性变化。

图 2.4：美国德克萨斯州费舍尔河的月度温度平均值和标准偏差

移动统计

在本节中，我们将继续研究 Fisher River 数据集。 通过在原始时间序列上滑动某个大小的窗口并汇总每个窗口的数据来计算移动或滚动统计信息。 这也称为在时间索引上卷积。 卷积操作的两个重要参数是窗口大小和步幅长度。 前者定义了作为聚合函数输入的时间单位的数量，而后者则定义了每次计算之间沿时间索引的间隔。 例如，假设k的窗口大小和长度l的步幅用于在时间序列 x [上计算函数f。 1] ，x [2] ，...，x [n] 具有N观察值。 在这种情况下，获得的运动统计量为 f（x [1] ，x [2] ，...，x [t）]，f（x [1 + 1] ，x [2 + 1] ，...，x [t + 1] ），依此类推。 请注意，每次通过将时间窗口向右移动l个时间单位来计算函数。

移动平均是函数f的特例，只需要简单地对时间窗口中的观测值求平均即可。

让我们演示如何在 Fisher River 数据集上计算移动平均值。 我们将计算每周移动平均值，该平均值将窗口大小设置为 7，并将窗口向右滑动一位：

weekly_moving_average = df['Mean temparature'].rolling(7).mean() 
Calculating monthly averages can be done as follows with a window of size thirty. 
monthly_moving_average = df['Mean temparature'].rolling(30).mean()

rolling函数仅将窗口大小作为参数。 因此，要增加一个以上的跨步长度，我们仍然如前所示计算移动平均值，然后对所得的序列进行切片以获得所需的输出。 对于跨越两个时间单位的步幅，我们使用以下代码：

weekly_moving_average_2stride = df['Mean temparature'].rolling(7).mean()[::2] 
monthly_moving_average_2stride = df['Mean temparature'].rolling(30).mean()[::2]

在时间序列分析中，最常基于步幅为 1 的步长移动统计信息，因此您几乎不需要滚动功能。 下图中绘制了原始数据以及每周和每月平均值，该图显示了移动平均值所产生的降噪效果和随后的平滑效果：

图 2.5：美国得克萨斯州费舍尔河的温度移动平均值

固定过程

诸如集中趋势，离散度，偏度和峰度之类的数据属性称为样本统计量。 均值和方差是最常用的两个样本统计量。 在任何分析中，数据都是通过从较大人群的样本中收集信息来收集的。 然后根据样本数据估算均值，方差和其他属性。 因此，这些被称为样本统计。

统计估计理论中的一个重要假设是，为了使样本统计可靠，总体不随样本中的个体或收集数据的时间发生任何根本或系统的变化。 此假设可确保样本统计信息不会改变，并且将适用于用于估计的样本之外的实体。

该假设也适用于时间序列分析，因此从简单估计中得出的均值，方差和自相关可以用作未来事件的合理估计。 在时间序列分析中，此假设称为平稳性，它要求序列的内部结构不随时间变化。 因此，平稳性要求均值，方差和自相关相对于实际观察时间是不变的。 理解平稳性的另一种方法是，该序列具有恒定的均值和恒定的方差，而没有任何可预测和重复的模式。

平稳时间序列的一个流行示例是零均值序列，它是从正态分布中均值为零的样本集合。 下图说明了零均值序列，该序列是从零均值和单位方差的正态分布采样的点生成的。 尽管从正态分布中依次采样点并将其绘制为时间序列，但各个观测值是独立的并且分布相同（iid）。 零均值序列没有显示任何时间模式，例如趋势，季节性和自相关。

图 2.6：零均值时间序列

但是，大多数实际时间序列不是固定的。 非平稳性主要是由于趋势和季节性的存在而影响不同时间点的均值，方差和自相关。 但是，值得注意的是，在定义平稳性时，我们并未包括周期性波动。 这是因为周期性变化引起的波峰和波谷不会以固定的间隔发生，因此只能用外生变量来解释。 通常，从长远来看，没有可预测模式的时间序列是固定的（不考虑外在因素作为解释变量！）。

时间序列分析中的关键步骤是通过特殊的数学运算来统计验证平稳性并使非平稳时间序列反平稳。 为此，我们将讨论用于检测平稳性的增强 Dickey-Fuller（ADF）测试，并描述用于使非平稳时间序列反平稳的微分方法。 差异可以消除趋势和季节性因素。 下一部分将讨论分解方法以开发复杂时间序列的趋势和季节性模型。

差异化

差分的基本思想是采用时间序列的连续出现之间的差异Δ x [t] = x [t] -x [t-1] 使得Δ x [t] 具有恒定的均值和方差，因此可以将其视为平稳序列。 ADF 测试基于区分原始时间序列的想法，因此，我们将在讨论各种类型的区分技术之后对其进行说明。

一阶微分

一阶差异意味着在时间序列的连续实现之间取差异，以使差异Δ x [t] 是不规则的变化，没有任何长期趋势或季节性。 上一章讨论的随机游动模型是后续随机变化的总和，由 x [t] = x [tl] +Є [t] ，其中Є [t] 是根据正态分布的零平均随机数。 随机游走的特点是顺序较长的向上或向下趋势。 此外，他们采取了无法预料的方向变化。 基于这些特征，随机行走是不固定的。 但是，随机游动的第一差（Δ x [t] 等于随机噪声Є [t]。因此残留残差 随机游动的一阶微分之后是零均值平稳序列，将一阶差分求出的变换时间序列表示为：

*x<sup class="calibre28">'</sup>[t] = x[t] - x[t-1]*

转换后的时间序列具有 N-1 个观测值，均值为零且方差恒定。 该模型假设起始值为 *x [1] =0。*起始值可以推广为 x [1] = 0。 不同的时间序列将是：

*x<sup class="calibre28">'</sup>[t] = x[t] - x[t-1] = c + Єt*

如果起始值 x [1] = 0 为正，则随机游走趋于向上漂移，而当其为负值时，则趋于向下漂移。

从前面的讨论中可以明显看出，一阶差分是独立的，并且具有恒定的均值和恒定的方差，因此没有自相关。

验证一阶微分是否已平稳化时间序列的快速方法是绘制 ACF 函数，并对差分序列运行 Ljung-Box 测试。 Ljung-Box 测试确定观察到的自相关是否具有统计显着性。 Ljung-Box 检验的零假设是时间序列由随机变化组成，并且缺乏可预测的自相关，而替代假设则表明观察到的自相关不是随机的。

让我们用时间序列道琼斯工业平均指数（DJIA）来说明 Ljung-Box 检验，该检验序列在上一章中用于说明自相关函数（ACF）。 我们首先获取 DJIA Close 值的一阶差，然后在下图中绘制得到的一阶差的时间序列：

图 2.7：道琼斯工业平均指数的时间序列及其一阶差异

接下来，针对不同的延迟计算 ACF，并通过 Ljung-Box 测试进行验证。 下图显示了 DJIA 收盘价的 ACF 以及一阶差的时间序列。 请注意，对于不同的序列，ACF 没有显示可预测的模式，并且突然下降至接近零。 此外，对于滞后= 10，检验的 p 值为 0.894，这使我们接受了 Ljung-Box 检验对于差分序列的零假设。

图 2.8：道琼斯工业平均指数及其一阶差异的自相关

现在，我们将简要讨论用于运行 DJIA 时间序列的 Ljung-Box 测试的代码。 假设原始时间序列已读入pandas.DataFrame djia_df，则一阶差异的计算如下：

first_order_diff = djia_df['Close'].diff(1)

对于 Ljung-Box 测试，我们使用statsmodels.tsa.stattools程序包中的acf功能。 acf函数用于返回自相关，置信区间，Q 统计量和测试的 p 值。

该示例的完整代码在本书的 GitHub 存储库的code文件夹下的 Jupyter 笔记本Chapter_2_First_Order_Differencing.ipynb中。

二阶微分

在某些情况下，一阶微分不会使时间序列平稳，因此数据会在另一个时间进行差分以生成固定的时间序列。 因此，二阶差分时间序列生成如下：

*x<sup class="calibre28">"</sup>[t] = x<sup class="calibre28">'</sup>[t] - x<sup class="calibre28">'</sup>[t-1] = (x[t] - x[t-1]) - (x[t-1] - x[t-2]) = x[t] - 2[xt-1]+x[t-2]*

由二阶微分产生的时间序列具有 N-2 个观测值。 几乎从来不需要执行高于二阶的阶数微分。

季节性差异

当一个时间序列具有已知时间段m时间索引的季节性时，可以通过取 x [t] 和 x [tm]。 这些滞后长度为m的差异表示一年中的季节或季度。 在这种情况下，m = 12，并且彼此之间相隔一年的原始观测值之间存在差异。 季节性差异可以表示为：

*x<sup class="calibre28">'</sup>[t] = x[t] - x[t-m] = Є[t]*

为了证明季节性差异的影响，我们将重新研究时间序列对费舍尔河日平均温度的影响。 我们已经看到了原始的时间序列和月平均值，两者显然都表现出强烈的季节性行为。 可以通过运行按月的分组操作（得出pandas.Series对象monthly_mean_temp）来获得月平均值。 使用pandas.plotting API 中的autocorrelation_plot函数计算和绘制该系列中的自相关，并在下图中显示。 autocorrelation_plot功能可用于检查时间序列中统计上显着的自相关的存在。

如下图所示，它还绘制了置信度分别为 95％（alpha = 0.05；粗虚线）和 99％（alpha = 0.01；细虚线）的上下置信区间。 费舍尔河月平均温度的 ACF 会在多个滞后的 99％置信区间上下波动。 。 费舍尔河月平均温度的 ACF 会在多个滞后的 99％置信区间上下波动。 因此，由于季节性因素，月平均温度形成了一个不稳定的时间序列。

图 2.9：美国得克萨斯州费舍尔河的月平均温度与时间序列的自相关

我们尝试通过采用以下季节差异来平稳每月平均值的时间序列：

seasonal_diff = monthly_mean_temp.diff(12)

代码的前一行返回seasonal_diff，它是pandas.Series。 季节性差异在其前 12 个元素中保留了空值，在进一步分析之前将其删除：

seasonal_diff = seasonal_diff[12:]

季节性差异似乎是随机变化，如下图所示：

图 2.10：美国德克萨斯州费舍尔河的月平均气温的季节性差异

我们再次使用autocorrelation_plot函数来生成差分序列的 ACF 和置信区间（置信度为 99％）。 在下图中我们可以看到，ACF 永远不会越过 99％的置信区间，而滞后从 0 到超过 35 不等：

图 2.11：美国德克萨斯州费舍尔河月平均温度的季节差异的自相关

通过对月平均数据运行stattools.acf函数，可以确认实际的 p 值，如下所示：

_, _, _, pval_monthly_mean = stattools.acf(monthly_mean_temp, unbiased=True, 
                                                                                        nlags=10, qstat=True, alpha=0.05) 
print('Null hypothesis is rejected for lags:',   np.where(pval_monthly_mean<=0.05))

前几行显示以下消息：

Null hypothesis is rejected for lags: (array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]),)

Ljung-Box 测试也在季节性差异序列上执行：

_, _, _, pval_seasonal_diff = stattools.acf(seasonal_diff, unbiased=True, 
                                                               nlags=10, qstat=True, alpha=0.05) 
print('Null hypothesis is rejected for lags:', np.where(pval_seasonal_diff<=0.05))

前面代码中的print函数返回以下行：

Null hypothesis is rejected for lags: (array([], dtype=int32),)

Ljung-Box 检验的零假设不会被滞后。

该示例的完整代码在本书的 GitHub 存储库的code文件夹下的 Jupyter 笔记本Chapter_2_Seasonal_Differencing.ipynb中。

在这一点上，重要的是要注意，在某些情况下，在运行季节性差异以实现转换后的数据的平稳性之后，会进行一阶差分。 所得的序列 x “ [t] 可以根据原始数据进行如下计算：

*x<sup class="calibre28">"</sup>[t] = x<sup class="calibre28">'</sup>[t] - x<sup class="calibre28">'</sup>[t-1] = (x[t] - x[t-m]) - (x[t-1] - x[t-m-1])*

可以通过探索性数据分析来确定差异策略的选择，就像到目前为止所描述的那样。 但是，当难以确定平稳所需的变换时，将执行 ADF 测试以提供明确的指导。

增强 Dickey-Fuller 测试

用于客观确定是​​否需要微分以使时间序列平稳的统计测试称为单位根测试。 我们讨论了几种此类测试，ADF 测试是单位根测试之一，最常用于验证原始时间序列中的非平稳性。 根据 ADF 测试，在存在自相关的情况下，原始序列的一阶差异 x ' [t] 可以表示为线性回归模型 时间索引和一阶差的最大差值，直到m时间索引的滞后。 x ' [t] 的线性回归公式如下：

在存在强自相关的情况下，原始序列需要区别，因此先前模型中的 Y x [t-1] 项在预测 时间索引 t-1 和t的变化。 因此，针对替代假设 H1 [：Y] < 0，ADF 的原假设为 H [0：Y] = 0。 换句话说，零假设是存在单位根或非平稳性，而另一种假设表明数据是平稳的。

对于给定的置信度，的检验统计量的原假设被拒绝。Y 小于临界负值。 如果序列已经是固定的，则线性回归模型表示一个随机漂移，其时间指数t的变化取决于先前的值，而不取决于先前的差，这对于平稳过程是同义的。

在假设很少需要高于三阶微分的微分来稳定时间序列的假设下，回归中要包括的滞后m的数目通常设置为 3。 确定m的一种实际方法是改变 mЄ{1,2，...，N'}，其中 *N'*是上限 限制滞后，并计算 ADF 检验的 p 值，直到在给定的置信度下拒绝零假设。

让我们应用 ADF 来验证 1963 年至 1970 年期间收集的美国航空每月飞机里程数据的平稳性。我们将使用statsmodels.tsa.stattools API 中的adfuller函数来运行测试。

在运行 ADF 测试之前，将时间序列加载到pandas.DataFrame中，并如下图所示进行绘制。 显然，时间序列具有上升趋势和季节性，因此是非平稳的，这将通过 ADF 测试进行验证。

图 2.12：1963-1969 年期间美国每月飞行航空里程的时间序列

读取和绘制时间序列的代码如下：

air_miles = pd.read_csv('datasets/us-airlines-monthly-aircraft-miles-flown.csv') 
air_miles['Air miles flown'].plot(ax=ax)

注意，我们已经使用了air_miles DataFrame的'Air miles flown'列的内置绘图功能。 现在，我们可以使用adfuller函数执行 ADF 测试，如下所示：

adf_result = stattools.adfuller(air_miles['Air miles flown'], autolag='AIC')

该函数的第一个参数是 DataFrame 的第二列，关键字参数autolag='AIC'通过最大化 Akaike 信息标准（AIC）。 或者，可以根据用户输入到关键字参数 maxlag 中的延迟数来运行测试。 我们更喜欢使用 AIC 而不是给出延迟，以避免在寻找运行测试所需的最佳延迟时避免反复试验。

毫不奇怪，ADF 测试的 p 值为 0.9945，这证实了我们对探索性数据分析的理解。 adfuller函数在一个元组中返回几个值。 我们已经使用adf_result变量存储结果，并引用其第二个元素adf_result[1]来检索测试的 p 值。 该函数返回的其他有趣的事物是usedlag，它是实际运行测试的滞后次数以及置信度为 1％，5％和 10％时测试统计量的临界值。

请参阅code/Chapter_2_Augmented_Dickey_Fuller_Test.ipynb以查看此示例的完整代码。

时间序列分解

时间序列分解的目的是对长期趋势和季节性建模，并结合它们来估计总体时间序列。 时间序列分解的两种流行模型是：

  • 加法模型
  • 乘法模型

加性模型将原始时间序列（ x [t] ）表示为趋势周期（ F [t] ）和季节性（ S [t] 的组成部分如下：

*x[t] = F[t] + S[t] + Є[t]*

调整趋势和季节分量后获得的残差Є [t] 是不规则变化。 当存在与时间相关的趋势周期成分，但独立的季节性不会随时间变化时，通常会应用加性模型。

当存在随季节变化的季节性时，将时间序列作为趋势，季节性和不规则成分的乘积的乘积分解模型非常有用：

*xt = F[t] x S[t] x Є[t]*

通过取对数，将乘法模型转换为各个分量的对数的加法模型。 乘法模型表示如下：

log（x [t] ）= log（F [t] ）+ log（S [t] ）+ log（Є [t] ）

在本节中，我们将讨论以下两种流行的方法来估计趋势和季节成分：

  • 均线法
  • 使用 Python 软件包statsmodels.tsa进行季节性和趋势分解

移动平均线

在本节中，我们将通过以下主题介绍移动平均线：

  • 移动平均值的计算及其在平滑时间序列中的应用
  • 使用移动平均线进行季节性调整
  • 加权移动平均线
  • 使用移动平均的时间序列分解

移动平均线及其平滑效果

在时间索引t上的移动平均值（MA）估计平均趋势周期成分 F [t]，并且是 通过对t±k 的时间段取平均值来计算，其中k是 MA 的范围：

采取移动平均值具有消除随机噪声来平滑原始时间序列的效果。 通常，观测的总数 m = 2k +1 用于将移动平均值描述为 m 阶 MA，此后将其表示为F_hat[t]^(m)。 让我们通过该示例来说明 1962 年至 1965 年 IBM 股票价格的移动平均线及其平滑效果。该时间序列的数据集可从这个页面下载。 原始时间序列（如下图蓝色所示）由于随机噪声而具有不规则运动。 以红色显示的 5 天（或 5 阶）MA 比原始序列更平滑，并且显示了趋势周期模式的估计值。 如下图所示，为期 5 天的移动平均显然对原始时间序列产生了一定的平滑作用。

该示例的代码在笔记本Chapter_2_Moving_Averages.ipynb的本书 GitHub 存储库的代码文件夹下。

原始时间序列已加载到panda.DataFrame中。 然后，为方便起见，将具有股票价格的第二列重命名为Close_Price。 接下来，通过将滚动函数应用于DataFrame的Close_Price列来计算 5 天 MA。 实现如下：

ibm_df = pd.read_csv('datasets/ibm-common-stock-closing-prices.csv') 
ibm_df.rename(columns={'IBM common stock closing prices': 'Close_Price'}, 
               inplace=True) 
ibm_df['5-Day Moving Avg'] = ibm_df['Close_Price'].rolling(5).mean() 
The plotting is done using the in-built plot function of pandas.Series objects. 
fig = plt.figure(figsize=(5.5, 5.5)) 
ax = fig.add_subplot(2,1,1) 
ibm_df['Close_Price'].plot(ax=ax) 
ax.set_title('IBM Common Stock Close Prices during 1962-1965') 
ax = fig.add_subplot(2,1,2) 
ibm_df['5-Day Moving Avg'].plot(ax=ax, color='r') 
ax.set_title('5-day Moving Average')

图 2.13：1962-1965 年期间 IBM 股票价格及其 5 天移动平均线

前面提到的具有奇数阶 m = 2 x 2 + 1 = 5 的 MA 对称，并且在时间索引t的两侧具有相等的观察数。 正在计算 MA 的时间。 然而，可能具有偶数阶 m = 2k 的不对称 MA。

偶数阶移动平均值的不对称性可以通过采用偶数阶的第二移动平均值来消除。 让我们通过考虑第一个移动平均数为 2 并计算如下来说明这一点：

另一个二阶移动平均线应用于序列F_hat[t]^(2)时，将产生最终对称移动平均线，该对称移动平均线在原始时间序列的索引t处的观测值的两侧都有一个观测值：

通常，我们可以通过首先采用阶数移动平均值，然后是N顺序移动平均值来创建n × F_hat[t]^(m)移动平均值。 然而，我们必须注意到，为了产生对称 MA，M和N应该是偶数或奇数。 我们将通过如下所示的3 × F_hat[t]^(3) MA 的计算来澄清这一点：

MA 可以按照原始序列进一步细分如下：

相反，让我们考虑可以用原始时间序列表示的2 × F_hat[t]^(3)，如下所示：

x [t] 在右侧的观察数量多于在左侧的观察数量，因此它是不对称的。 顺便说一句，对称 MA 在估计季节性时间序列的趋势周期方面具有有趣的应用。 下一部分将对此进行说明。

还值得注意的是，形式为n × F_hat[t]^(m)的重复移动平均值产生了来自原始时间序列的观测值的加权平均值。 在时间索引t上分配给观察的权重最高，并且随着我们从t的两侧进一步移开，权重下降。 因此，n × F_hat[t]^(m)可用于生成加权移动平均值。

此处给出的图显示了 IBM 股价在最初 45 天之内的六个移动平均线，即F_hat[t]^(2)，2 × F_hat[t]^(2)，F_hat[t]^(4)，2 × F_hat[t]^(4)，F_hat[t]^(3)和3 × F_hat[t]^(3)。 如图所示，得到的 MA 的平滑度随着以m和重复次数n的顺序增加而提高。

在下图中，以实线给出了F_hat[t]^(m)类型的 MA，而以虚线绘制了n × F_hat[t]^(m)类型的 MA。

为了生成上述六个 MA，我们广泛使用了rolling和mean函数，如以下代码段所示：

MA2 = ibm_df['Close_Price'].rolling(window=2).mean() 
TwoXMA2 = MA2.rolling(window=2).mean() 

MA4 = ibm_df['Close_Price'].rolling(window=4).mean() 
TwoXMA4 = MA4.rolling(window=2).mean() 

MA3 = ibm_df['Close_Price'].rolling(window=3).mean() 
ThreeXMA3 = MA3.rolling(window=3).mean()

图 2.14：IBM 股票价格的不同移动平均线

使用移动平均线进行季节性调整

n × F_hat[t]^(m)移动平均值的加权平均属性已应用于具有季节性的平滑数据，以便生成趋势周期的估计值。 例如，给定季度观测值，我们可以应用2 × F_hat[t]^(4) MA 来平滑数据。 要了解其工作原理，让我们根据原始时间序列展开2 × F_hat[t]^(4)：

对于季度数据，2 × F_hat[t]^(4) MA 的第一项和最后一项与同一季度相对应，但是在连续的几年中，并且取其季度平均值，该变化在一年中会有所不同。 同样，对于月度数据，2 × F_hat[t]^(12)将生成一个平滑序列作为趋势周期分量的估计。 对于每周观测等具有周期性的数据，以F_hat[t]^(2k + 1) MA（例如F_hat[t]^(7)）作为每周数据可消除季节性。

使用来自澳大利亚的啤酒生产数据说明了使用2 × F_hat[t]^(2k)形式的加权 MA 的方法。 该时间序列表示从 1956 年 3 月到 1994 年 6 月的每个季度的啤酒产量。以下是原始时间序列和一系列2 × F_hat[t]^(4) MA：

图 2.15：1957-1992 年期间的季度啤酒产量及其 2 比 4 的移动平均值

读取数据集并计算移动平均值的代码如下：

beer_df = pd.read_csv('datasets/quarterly-beer-production-in-aus-March 1956-     June 1994.csv') 
beer_df.rename(columns={'Quarterly beer production in Australia: megalitres.  March 1956 ? June 1994': 'Beer_Prod'}, 
               inplace=True 
              ) 
MA4 = beer_df['Beer_Prod'].rolling(window=4).mean() 
TwoXMA4 = MA4.rolling(window=2).mean() 
TwoXMA4 = TwoXMA4.loc[~pd.isnull(TwoXMA4)]

季度啤酒产量的原始时间序列具有趋势和季节性，因此不是固定的。 让我们看看是否可以通过先删除趋势分量然后采用季节性差异来平稳时间序列。

我们首先从去除趋势分量后剩余的残差开始：

residuals = beer_df['Beer_Prod']-TwoXMA4 
residuals = residuals.loc[~pd.isnull(residuals)]

下图显示了删除趋势周期分量后剩余的residuals：

图 2.16：减去季度啤酒产量与其移动平均值之差后的剩余残留量

在这一点上，我们将通过绘制自相关函数以及 99％的置信区间来检查残差是否已经平稳（虽然不太可能！）。 要获得此图，如此处所示，我们将使用pandas.plotting API 中的autocorrelation_plot函数：

autocorrelation_plot(residuals, ax=ax)

图 2.17：季度啤酒产量的自相关

显然，对于几个滞后值，残差与 ACF 跳到置信区间之外具有很强的自相关性。 因此，我们需要对残差采取季节性差异。 可以根据以下事实确定季节性周期：原始数据是从年份的所有季度获得的，并显示该季度的季节性。 这意味着一年第一季度的残差在大小上与前一年及后续年份的第一季度的残差接近。 该观察结果使我们在四个时间单位的时间段内得出的差异如下：

residuals_qtr_diff = residuals.diff(4) 
residuals_qtr_diff = residuals_qtr_diff.loc[~pd.isnull(residuals_qtr_diff)]

我们期望residuals_qtr_diff是一系列随机变化，没有季节性和可预测的自相关。 为了验证是否是这种情况，我们在residuals_qtr_diff上运行autocorrelation_plot函数，并获得以下具有随机 ACF 的图。 此外，ACF 仅两次落后就落在 99％置信区间之外。 这意味着，通过季节性差异使残差平稳。

图 2.18：减去季度啤酒产量与其移动平均值之差后剩余残留物的自相关

我们固定季度啤酒生产时间序列的方法可以总结如下：

  1. Take seasonal 2 × F_hat[t]^(4) MA.

  2. 通过从原始时间序列中删除2 × F_hat[t]^(4) MA 来计算残差。

  3. 检查残差的 ACF 的随机性。

  4. 如果残差的 ACF 已经是随机的，则残差是固定的，如果没有，则继续下一步。

  5. 取周期性为 4 的残差的季节性差异，并检查差异序列的 ACF 的随机性。

因此，季节性移动平均数和残差的季节性差异基本上使原始时间序列平稳。

加权移动平均线

在前面的部分中，我们将n × F_hat[t]^(m)形式的移动平均值表示为来自原始时间序列的观测值的加权和。 我们看到，观测值的权重是如何下降的，而不是与正在计算 MA 的时间索引相距甚远。 加权移动平均值的概念是对称的，可以将其推广到任何时间序列应用，如下所示：

权重w[t-k] + ... + w[t+k] = 1。

对于简单的F_hat[t]^(m)，权重为1/m。 已经观察到，加权 MA 可以更平滑地估计趋势周期。 这归因于以下事实：在原始时间序列中用于计算加权 MA 的所有观测值都没有像n × F_hat[t]^(m) MA 系列那样被赋予相等的权重。 接近时间索引t的观测值被分配较高的权重，而距离较远的观测值的权重较小。 这种根据与感兴趣的时间索引的接近程度分配不同权重的方法可以更好地平滑数据。

使用移动平均的时间序列分解

在分析季度啤酒生产数据时，我们开发了一种使用季节性平均价格和季节性差异对非平稳时间序列进行平稳的方法。 我们使用季节性移动平均线作为趋势周期成分的估计值，并计算出 MA 留下的残差的周期性差异。

另一种方法可能是从原始序列中扣除季节均值和季节残差，并检查最终残差的随机性。 该方法假定啤酒生产系列是趋势周期和季节成分的加总和，并且除去上述两个变量后剩下的是随机变化。 实际上，可以用移动平均线以这种方式分解时间序列。 让我们在本节中详细说明移动平均驱动时间序列分解。

形式为n ×F_hat[t]^(m)的移动平均值具有平滑原始时间序列并估计趋势周期分量的特性。m和n的选择对于确定趋势周期成分至关重要。 通常，m是季节性数据的周期性，并且是先验的或通过探索性数据分析确定。 如果m是偶数，则F_hat[t]^(m)移动平均值将是不对称的，因此在连续出现的季节中缺少必要的平滑特性。 因此，第二次采用移动平均线来生成2 × F_hat[t]^(m)序列，该序列具有必要的季节性平滑度，并且可以更好地估计趋势周期分量。 但是，在奇数周期的情况下，F_hat[t]^(m)被用作趋势周期分量的估计。 因此，取决于 m 是偶数还是奇数，趋势周期分量的估计是F_hat[t] = 2 × F_hat[t]^(m)或F_hat[t]^(m)。

在使用移动平均值的时间序列分解中，假定季节分量在一年中或一周或一周中是恒定的。 通过对趋势周期进行调整后剩余的残值按季节进行平均，可以得出季节性成分的估计值。 例如，如果数据是每月数据，则我们将趋势周期调整的残差按月取平均值。

加法模型和乘法模型之间的计算会略有不同。 例如，在通过应用适当的移动平均线估算趋势周期成分F_hat[t]时，x[t] - F_hat[t]计算出加性模型的残差，而对于乘法模型，其残差为x[t]/F_hat[t]。 现在，从残差中估算季节成分S_hat[t]，作为季节平均值。

最后，通过对原始序列进行趋势周期和季节性调整，可以得出不规则的变化，如下所示：

  • 对于加法分解模型：ε[t] = x[t] - F_hat[t] - S_hat[t]
  • 对于乘法分解模型：ε[t] = x[t] / (F_hat[t] × S_hat[t])

让我们用美国航空每月飞行的密尔飞机按时间序列的移动平均值来说明这种时间序列分解的方法。 我们对此数据集应用加法和乘法模型。

我们首先将数据集读取到pandas.DataFrame air_miles中。 接下来，通过2 × F_hat[t]^(12)移动平均值估算趋势周期成分，如下所示：

MA12 = air_miles['Air miles flown'].rolling(window=12).mean() 
trendComp = MA12.rolling(window=2).mean()

对于加性模型，通过从原始时间序列中减去趋势周期并获取残差的按月平均值来获得季节性成分。 分组操作用于计算季节性分量，如下所示：

residuals = air_miles['Air miles flown'] - trendComp 

#To find the sesonal compute we have to take monthwise average of these residuals 
month = air_miles['Month'].map(lambda d: d[-2:]) 
monthwise_avg = residuals.groupby(by=month).aggregate(['mean']) 
#Number of years for which we have the data 
nb_years = 1970-1963+1 

seasonalComp = np.array([monthwise_avg.as_matrix()]*nb_years).reshape((12*nb_years,))

给定趋势周期和季节成分的估计值，对于加性模型，我们得到如下不规则变化：

irr_var = air_miles['Air miles flown'] - trendComp - seasonalComp

下图显示了原始时间序列，趋势周期，季节性和不规则部分：

图 2.19：每月航空里程上的时间序列的加法分解

在irr_var上进行的 ADF 测试得出 p 值为 0.0658。 在 90％的置信度下（alpha = 0.10），可以接受关于不规则变化的平稳性的零假设。 但是，让我们尝试通过乘法模型进行进一步的改进。

MA 的趋势周期估计甚至适用于乘法模型。 但是，季节性成分的计算会发生变化：

residuals = air_miles['Air miles flown'] / trendComp

请注意，在乘法模型的情况下，季节性residuals从原始时间序列划分趋势周期组件。 但是，通过对残差的逐个操作保持不变。

最后，我们得出乘法模型的不规则变异，如下所示：

irr_var = air_miles['Air miles flown'] / (trendComp * seasonalComp)

从乘法模型获得的不规则变化的 ADF 测试得出的 p 值约为 0.00018，这比从加法模型获得的 p 值小得多。

鼓励读者在 Jupyter 笔记本code/Chapter_2_Time_Series_Decomposition_by_Moving_Averages.ipynb中查看本节中完成的分析的完整代码。

下图绘制了原始时间序列以及趋势周期，季节性和不规则成分。 注意，季节和不规则分量的数量级要比从加法模型获得的分量小：

图 2.20：每月航空里程的时间序列的乘法分解

使用 statsmodels.tsa 进行时间序列分解

到目前为止，我们已经讨论了如何使用 MA 来估计时间序列的趋势周期和季节成分。 MA 的方法是在一个简单的假设下工作的，即季节性变化在连续的几年，几周或一段适合给定用例的时期内是恒定的。 但是，对于需要高级方法的一些应用来说，恒定的季节性可能是有效的，例如使用“局部加权平滑散点图”的季节性和趋势分解（通常也称为 STL 方法）。

在本节中，我们将使用 Python statsmodels.tsa包处理具有复杂模式的时间序列。 我们的目标是估计趋势周期和季节性因素。 此外，我们将使用趋势周期和季节估计来平稳化将由 ADF 测试验证的数据。

让我们考虑一下美国威斯康星州每月就业的时间序列，以说明上述方法。 数据为 1961 年 1 月至 1975 年 10 月的数据，可从这个页面下载。

我们首先将数据集加载到pandas.DataFrame wisc_emp中，然后对原始时间序列运行 ADF 测试，如下所示：

wisc_emp = pd.read_csv('datasets/wisconsin-employment-time-series.csv') 
adf_result = stattools.adfuller(wisc_emp['Employment'], autolag='AIC')

ADF 测试在每月就业序列上的高 p 值为 0.9810，表明原始时间序列是不稳定的。 因此，我们尝试分解时间序列，然后使用statsmodels.tsa API 中的seasonal.seasonal_decompose函数使其平稳。 让我们首先尝试加法模型进行分解：

decompose_model = seasonal.seasonal_decompose(wisc_emp.Employment.tolist(), freq=12, model='additive')

seasonal.seasonal_decompose中的参数freq是季节性行为的周期性，而原始时间序列是每月观测值，我们怀疑该周期性为 12，可以通过探索性数据分析来验证。

可通过seasonal.seasonal_decompose函数返回的对象decompose_model的属性访问分解时间序列的趋势周期，季节和残差分量。 这些组件可以从 decompose_model 的以下属性中获取：

  • decompose_model.trend-趋势周期组件
  • decompose_model.seasonal-季节性成分
  • decompose_model.resid-不规则的变化

现在，我们对加性模型的残差进行 ADF 测试，并获得 0.00656 的 p 值，该值比从原始时间序列获得的 p 值小得多。 但是，我们还将建立一个乘法模型：

decompose_model = seasonal.seasonal_decompose(wisc_emp.Employment.tolist(), freq=12, model='multiplicative')

ADF 对乘法分解残差的测试得出的 p 值为 0.00123，甚至比加法分解得到的 p 值还要小。 在置信水平为 99％（alpha = 0.01）的情况下，我们可以拒绝 ADF 检验的原假设，并得出结论，乘法分解模型的残差是平稳序列。

运行本节中示例的代码在 Jupyter 笔记本code/Chapter_2_Time_Series_Decomposition_using_statsmodels.ipynb中。

下图中绘制了原始时间序列以及各个组成部分：

图 2.21：使用 statsmodels.tsa API 分解每月就业的时间序列

概括

在本章的开头，我们讨论了高级数据处理技术，例如重采样，分组依据和移动窗口计算，以从时间序列中获取汇总统计信息。 接下来，我们描述了平稳的时间序列，并讨论了假设的统计检验，例如 Ljung-Box 检验和增强 Dickey Fuller 检验，以验证时间序列的平稳性。 平稳化非平稳时间序列对于时间序列预测很重要。 因此，我们讨论了两种平稳时间序列的方法。 首先，已经描述了覆盖第一，第二和季节差分的差分方法，用于使非平稳时间序列平稳。 其次，讨论了使用statsmodels.tsa API 进行加性和乘性模型的时间序列分解。 在下一章中，我们将深入研究处理噪声时间序列数据的指数平滑技术。