时间序列分析实战-全-

时间序列分析实战（全）

原文：Practical Time Series Analysis

协议：CC BY-NC-SA 4.0

零、前言

本书是有关使用 Python 进行时间序列分析的简介。我们旨在为您提供该学科基本概念的清晰概述，并描述适用于该行业中常见的分析用例的有用技术。由于有太多的项目需要根据过去的数据进行趋势分析和预测，因此时间序列分析是任何现代数据科学家的知识宝库中的重要工具。本书将为您提供工具和技术，使您可以自信地思考问题并提出时间序列预测中的解决方案。

为什么选择 Python？ Python 正在迅速成为跨不同行业的数据科学项目的首选。大多数最先进的机器学习和深度学习库都具有 Python API。结果，许多数据科学家更喜欢使用 Python 来实现整个项目管道，该管道包括数据整理，模型构建和模型验证。此外，Python 提供了易于使用的 API 来处理，建模和可视化时间序列数据。此外，Python 是开发 Web 应用后端的一种流行语言，因此吸引了广大软件专业人员。

现在，让我们看看您可以从本书的每一章中学到什么。

本书涵盖的内容

第 1 章和“时间序列简介”首先讨论三种不同类型的数据集-横截面，时间序列和面板。讨论了从横截面到时间序列的过渡以及数据分析增加的复杂性。描述了使时间序列数据变得特别的特殊数学属性。几个示例演示了如何使用探索性数据分析来可视化这些属性。

第 2 章，“了解时间序列数据”涉及三个主题，即通过重采样，分组依据对时间序列数据进行高级预处理和可视化 ]和移动平均线的计算；平稳性和统计假设检验，以检测时间序列中的平稳性；以及用于平稳化非平稳时间序列的各种时间序列分解方法。

第 3 章，“基于指数平滑的方法”涵盖了使用 Holt-Winters 方法进行一阶捕获级别，二阶平滑级别和趋势以及高阶平滑的基于平滑的模型。说明了该图，该图捕获了时间序列数据集中的水平，趋势和季节性。

第 4 章和“自回归模型”讨论了用于预测的自回归模型。本章涵盖移动平均值（MA），自回归（AR），自回归移动平均值（ARMA）和自回归综合移动平均值（ARIMA）的详细实现，以捕获预测期间时间序列数据中不同级别的干扰。

第 5 章，“用于时间序列预测的深度学习”讨论了可以直接用于开发时间序列数据的预测模型的最新深度学习算法。递归神经网络（RNN）是对数据序列进行建模的自然选择。在本章中，将介绍不同的 RNN，例如 Vanilla RNN，门控循环单元和长期短期记忆单元，以开发时间序列数据的预测模型。在概念上讨论了开发这些 RNN 的数学公式。案例研究使用 Python 的“ keras”深度学习库解决。

附录和 Python 入门，您将快速，轻松地了解 Python。如果您不熟悉 Python 或正在寻找如何开始使用编程语言，阅读本附录将帮助您克服最初的障碍。

这本书需要什么

您将需要 Anaconda Python 发行版来运行本书中的示例，并编写自己的 Python 程序以进行时间序列分析。可从这个页面免费下载。

本书的代码示例是使用 Jupyter Notebook 开发环境编写的。要运行 Jupyter Notebook，您需要安装 Anaconda Python Distribution，它具有 Python 语言要点，解释器，用于开发示例的包以及 Jupyter Notebook 服务器。

这本书是给谁的

本书的主题预期对以下人员有用：

数据科学家，具有统计，机器学习以及模型构建和验证背景的专业人员
数据工程师，具有软件开发背景的专业人员
希望在生成数据驱动的业务见解方面发展专业知识的软件专业人员

约定

在本书中，您将找到许多可以区分不同类型信息的文本样式。以下是这些样式的一些示例，并对其含义进行了解释。

代码块设置如下：

import os
import pandas as pd
%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns

文本代码以字体和颜色突出显示，如下所示：pandas.DataFrame。文件和文件夹的名称也以相同的样式显示，例如Chapter_1_Models_for_Time_Series_Analysis.ipynb or datasets/DJIA_Jan2016_Dec2016.xlsx

在本书的多个地方，我们都引用了外部 URL 来引用数据集或其他信息的来源。 URL 将以以下文本样式显示

新术语和重要词以粗体显示。您在屏幕上看到的字词，例如在菜单或对话框中，将以如下形式显示在文本中：“为了下载新模块，我们将转到“文件” |“设置” |“项目名称” |“项目解释器”。”

警告或重要提示如下所示。

提示和技巧如下所示。

读者反馈

始终欢迎读者的反馈。让我们知道您对这本书的看法-您喜欢或不喜欢的书。读者反馈对我们很重要，因为它可以帮助我们开发出您真正能充分利用的标题。要向我们发送一般性反馈，只需通过电子邮件发送feedback@packtpub.com，然后在您的消息主题中提及该书的标题。如果您有专业知识的主题，并且对写作或撰写书籍感兴趣，请参阅这个页面上的作者指南。

客户支持

既然您是 Packt 书的骄傲拥有者，我们可以通过很多方法来帮助您从购买中获得最大收益。

下载示例代码

您可以从这个页面的帐户中下载本书的示例代码文件。如果您在其他地方购买了此书，则可以访问这个页面并注册以将文件直接通过电子邮件发送给您。您可以按照以下步骤下载代码文件：

使用您的电子邮件地址和密码登录或注册到我们的网站。
将鼠标指针悬停在顶部的“支持”选项卡上。
单击“代码下载和勘误”。
在搜索框中输入书籍的名称。
选择您要下载其代码文件的书。
从购买本书的下拉菜单中选择。
单击代码下载。

下载文件后，请确保使用以下最新版本解压缩或解压缩文件夹：

Windows 的 WinRAR / 7-Zip
适用于 Mac 的 Zipeg / iZip / UnRarX
适用于 Linux 的 7-Zip / PeaZip

本书的代码包也托管在 GitHub 上。我们还从这个页面提供了丰富的书籍和视频目录中的其他代码包。去看一下！

勘误

尽管我们已尽一切努力确保内容的准确性，但还是会发生错误。如果您发现我们的其中一本书中有错误-可能是文字或代码中的错误-请向我们报告，我们将不胜感激。这样，您可以使其他读者免于沮丧，并帮助我们改进本书的后续版本。如果您发现任何勘误，请访问这个页面进行报告，选择您的图书，点击“勘误提交表格”链接，然后输入勘误的详细信息。勘误一经确认，您的提交将被接受，勘误将被上载到我们的网站或添加到该标题的勘误部分下的任何现有勘误列表中。要查看以前提交的勘误，请转到这个页面，然后在搜索字段中输入书籍的名称。所需信息将出现在“勘误”部分下。

海盗行为

互联网上版权材料的盗版在所有媒体中都是一个持续存在的问题。在 Packt，我们非常重视版权和许可的保护。如果您在互联网上以任何形式发现我们的作品的任何非法副本，请立即向我们提供位置地址或网站名称，以便我们寻求补救。请通过copyright@packtpub.com与我们联系，并提供指向可疑盗版材料的链接。感谢您在保护我们的作者方面的帮助以及我们为您带来有价值的内容的能力。

问题

如果您对本书的任何方面都有疑问，可以通过questions@packtpub.com与我们联系，我们将尽最大努力解决该问题。

一、时间序列简介

近年来，目睹了统计和机器学习的广泛应用，以从几乎所有工业领域的数据中得出可行的见解和商业价值。因此，业务分析师和软件专业人员必须能够处理不同类型的数据集。通常，数据是时间序列，以一系列有关系统或过程的定量观察的形式出现，并在连续的时间点进行。通常，时间点之间的距离是相等的。时间序列数据的示例包括国内生产总值，销售量，股票价格，在数年，数月，数天，数小时等的时间范围内记录的天气属性。观察的频率取决于变量的性质及其应用。例如，每年公开报告用于衡量一个国家的年度经济进步的国内生产总值。销量每月，每季度或每两年发布一次，尽管较长时间的数字可能是通过汇总更细粒度的数据（例如每日或每周的销量）而生成的。每秒都有关于股票价格和天气属性的信息。另一方面，有几个物理过程会在几分之一秒内生成时间序列数据。

成功利用时间序列数据将导致随着时间的推移监视系统的运行状况。例如，从其季度利润率跟踪公司的业绩。时间序列分析旨在将这些数据用于多种目的，可以大致归类为：

理解和解释随着时间推移产生观察到的系统或过程状态的潜在力量
根据可观察的特征预测系统或过程的未来状态

为了实现上述目标，时间序列分析应用了不同的统计方法来探索和建模时间序列数据的内部结构，例如趋势，季节性波动，周期性行为和不规则变化。存在几种数学技术和编程工具可以有效地设计计算机程序，这些程序可以探索，可视化和建模时间序列数据中的模式。

但是，在深入研究这些技术之前，本章旨在解释以下两个方面：

时间序列数据与非时间序列数据之间的差异
时间序列的内部结构（在上一段中已经简要提及了其中的一些）

对于解决问题，读者会发现本章有用以：

区分时间序列数据和非时间序列数据，因此选择正确的方法来表述和解决给定的问题。
为时间序列问题选择适当的技术。根据应用的不同，可以选择将重点放在时间序列数据的一个或多个内部结构上。

在本章的最后，您将了解您在分析项目中可能必须处理的不同类型的数据集，并能够区分时间序列与非时间序列。您还将了解使数据成为时间序列的特殊内部数据结构。从本章中学到的总体概念将有助于选择处理时间序列的正确方法。

本章将涵盖以下几点：

了解您在分析项目中可能遇到的不同类型的数据
了解构成时间序列的数据的内部结构
自相关是时间序列的最重要的内部结构，通常是时间序列分析的主要重点

不同类型的数据

业务分析师和数据科学家在其分析项目中会遇到许多不同类型的数据。在学术和工业项目中通常可以找到的大多数数据可以大致分为以下几类：

横截面数据
时间序列数据
面板数据

了解解决问题所需的数据类型以及可以从可用资源中获取哪种类型的数据对于提出问题和选择正确的分析方法非常重要。

横截面数据

人口的横截面数据或横截面是通过在同一时间点从多个人的观察中获得的。横截面数据可以包含在不同时间点获得的观察结果，但是，在这种情况下，时间本身在分析中没有任何重要作用。特定年份中高中生的 SAT 分数是横断面数据的一个示例。给定年份中国家的国内生产总值是横截面数据的另一个示例。用于客户流失分析的数据是横截面数据的另一个示例。请注意，在学生的 SAT 分数和国家/地区 GDP 的情况下，所有观测值都是在一年内进行的，这使两个数据集具有横截面。本质上，在两种情况下，横截面数据都代表给定时间实例的快照。但是，可以从诸如数年和数月之类的时间范围内获得用于客户流失分析的客户数据。但是出于分析目的，时间可能不会发挥重要作用，因此尽管客户流失数据可能来自多个时间点，但仍可以将其视为横截面数据集。

通常，对横截面数据的分析从绘制变量图开始，以可视化它们的统计属性，例如集中趋势，离散度，偏度和峰度。下图以 2010 年军费开支占 85 个国家的国内生产总值百分比的单变量示例对此进行了说明。通过采用一年的数据，我们可以确保其横截面性质。该图结合了标准化的直方图和核密度图，以突出显示军事费用数据的不同统计属性。

从该图可以明显看出，军事支出略有偏斜，大约在 1.0％左右达到最高峰。在 6.0％和 8.0％附近也可以观察到几个次要峰。

图 1.1：单变量横截面数据示例

还可以对多个变量进行探索性数据分析（例如上图中的数据分析），以了解它们的联合分布。让我们通过考虑这些国家的中央政府的总债务及其 2010 年的军费开支来说明双变量分析。下图以核密度图的形式显示了这些变量的联合分布。双变量联合分布表明两者之间没有明显的相关性，只是军事支出和中央政府债务的价值较低。

图 1.2：双变量横截面数据示例

值得注意的是，如前面的示例所示，横截面数据的分析超出了探索性数据分析和可视化的范围。诸如横截面回归之类的高级方法在几个解释变量和因变量之间拟合了线性回归模型。例如，在客户流失分析的情况下，目标可能是在由搅动或不搅动描述的顾客属性和顾客行为之间拟合逻辑回归模型。逻辑回归模型是离散和二进制结果的广义线性回归的一种特殊情况。它解释了导致客户流失的因素，并可以预测新客户的结果。由于时间不是此类横截面数据中的关键因素，因此可以在将来的某个时间点为新客户获得预测。在本书中，我们讨论了对时间序列数据进行建模的技术，其中时间和观测的顺序性质是进行分析的关键因素。

有关国家军事支出和国债的示例数据集已从世界银行的开放数据目录中下载。您可以在本书的 GitHub 存储库的datasets文件夹下的WDIData.csv文件中找到数据。

本书中的所有示例均附带有相同的 Python 实现。因此，现在让我们讨论为生成上述图形而编写的 Python 程序。在绘制图形之前，我们必须将数据集读入 Python，并根据数据集中发现的列和行熟悉数据的基本结构。本书中用于示例和图形的数据集为 Excel 或 CSV 格式。我们将使用pandas包来读取和处理数据。为了可视化，使用了matplotlib和seaborn。让我们首先导入所有程序包以运行此示例：

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

print_function已从__future__包中导入，以使print可以作为可能使用 Python 2.x 版本的读者的功能。在 Python 3.x 中，默认情况下print是一个函数。由于此代码是从 IPython 笔记本编写和执行的，因此%matplotlib inline可确保正确导入图形包并使其可在笔记本的 HTML 环境中工作。 os程序包用于如下设置工作目录：

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

现在，我们从 CSV 文件中读取数据并显示有关其的基本信息：

data = pd.read_csv('datasets/WDIData.csv') 
print('Column names:', data.columns)

这为我们提供了以下输出，显示了数据集的列名：

Column names: Index([u'Country Name', u'Country Code', u'Indicator Name', 
     u'Indicator Code', u'1960', u'1961', u'1962', u'1963', u'1964',  u'1965', 
       u'1966', u'1967', u'1968', u'1969', u'1970', u'1971', u'1972', u'1973', 
       u'1974', u'1975', u'1976', u'1977', u'1978', u'1979', u'1980', u'1981', 
       u'1982', u'1983', u'1984', u'1985', u'1986', u'1987', u'1988', u'1989', 
       u'1990', u'1991', u'1992', u'1993', u'1994', u'1995', u'1996', u'1997', 
       u'1998', u'1999', u'2000', u'2001', u'2002', u'2003', u'2004', u'2005', 
       u'2006', u'2007', u'2008', u'2009', u'2010', u'2011', u'2012', u'2013', 
       u'2014', u'2015', u'2016'], 
      dtype='object')

通过运行以下行，我们还可以从行和列数的角度了解数据的大小：

print('No. of rows, columns:', data.shape)

这将返回以下输出：

No. of rows, columns: (397056, 62)

该数据集具有近 40 万行，因为它捕获了 264 个不同国家的 1504 个世界发展指标。通过运行以下四行，可以获取有关指标和国家/地区的唯一数目的信息：

nb_countries = data['Country Code'].unique().shape[0] 
print('Unique number of countries:', nb_countries)

从数据结构上看，每一行都给出了有关由Indicator Name和Indicator Code列标识的指标以及对于国家/地区的观察值，该指标由Country Name和Country Code列表明。列1960至2016具有相同时间段内的指标值。基于对DataFrame中数据布局方式的了解，我们现在可以提取与可视化相关的行和列。

让我们首先准备另外两个DataFrames，它们得到与所有国家的指标Total Central Government Debt (as % of GDP)和Military expenditure (% of GDP)相对应的行。这是通过将原始DataFrame切片如下来完成的：

central_govt_debt = data.ix[data['Indicator Name']=='Central government debt, total (% of GDP)'] 
military_exp = data.ix[data['Indicator Name']=='Military expenditure (% of GDP)']

前两行创建两个新的DataFrames，即central_govt_debt和military_exp。通过运行以下两行，可以快速检查这些DataFrames的形状：

print('Shape of central_govt_debt:', central_govt_debt.shape) 
print('Shape of military_exp:', military_exp.shape)

这些行返回以下输出：

Shape of central_govt_debt: (264, 62) 
Shape of military_exp: (264, 62)

这些DataFrames具有我们需要的所有信息。为了在上图中绘制单变量和双变量横截面数据，我们需要2010列。在实际运行绘图代码之前，让我们快速检查2010列是否丢失。这是通过以下两行完成的：

central_govt_debt['2010'].describe() 
military_exp['2010'].describe()

分别产生以下输出：

count     93.000000 
mean      52.894412 
std       30.866372 
min        0.519274 
25%             NaN 
50%             NaN 
75%             NaN 
max      168.474953 
Name: 2010, dtype: float64 
count    194.000000 
mean       1.958123 
std        1.370594 
min        0.000000 
25%             NaN 
50%             NaN 
75%             NaN 
max        8.588373 
Name: 2010, dtype: float64

这就告诉我们 describe 函数无法计算第 25 个^{，第^个和第 75 ^个的四分位数，因此有一些遗漏的值需要避免。}

另外，我们希望Country Code列成为行索引。因此，执行以下几行：

central_govt_debt.index = central_govt_debt['Country Code'] 
military_exp.index = military_exp['Country Code']

接下来，我们通过从central_govt_debt和military_exp中获取非空的2010列来创建两个pandas.Series。然后将新创建的Series对象合并为一个DataFrame：

central_govt_debt_2010 = central_govt_debt['2010'].ix[~pd.isnull(central_govt_debt['2010'])] 
military_exp_2010 = military_exp['2010'].ix[~pd.isnull(military_exp['2010'])] 
data_to_plot = pd.concat((central_govt_debt_2010, military_exp_2010), axis=1) 
data_to_plot.columns = ['central_govt_debt', 'military_exp'] 
data_to_plot.head()

前几行返回下表，该表表明并非所有国家都掌握 2010 年中央政府债务和军事费用的信息：

| | central_govt_debt | Military_exp |
| 美国空军 | NaN | 1.897473 |
| 前 | NaN | 4.244884 |
| 白色的 | NaN | 1.558592 |
| ARB | NaN | 5.122879 |
| 是 | NaN | 6.119468 |
| 生气的 | NaN | 0.814878 |
| 手臂 | NaN | 4.265646 |
| ATG | 75.289093 | NaN |
| 出去 | 29.356946 | 1.951809 |
| 或者 | 79.408304 | 0.824770 |

要作图，我们只需要选择既有中央政府债务又有军事费用的国家。运行以下行，以筛选出缺少值的行：

data_to_plot = data_to_plot.ix[(~pd.isnull(data_to_plot.central_govt_debt)) & (~pd.isnull(data_to_plot.military_exp)), :]

通过运行以下行，可以显示已过滤的DataFrame的前五行：

data_to_plot.head()

	central_govt_debt	Military_exp
出去	29.356946	1.951809
或者	79.408304	0.824770
AZE	6.385576	2.791004
贝尔	7.022605	1.084631
BGR	21.286254	1.765384
出去	29.356946	1.951809
或者	79.408304	0.824770
AZE	6.385576	2.791004
贝尔	7.022605	1.084631
BGR	21.286254	1.765384

上表只有非空值，我们现在可以为横截面数据生成图了。以下代码行生成有关军事费用的单变量横截面数据的图：

plt.figure(figsize=(5.5, 5.5)) 
g = sns.distplot(np.array(data_to_plot.military_exp), norm_hist=False) 
g.set_title('Military expenditure (% of GDP) of 85 countries in 2010')

该图将另存为[GitHub]存储库的plots/ch1文件夹下的png文件。我们还将通过运行以下代码来生成军费与中央政府债务之间的二元图：

plt.figure(figsize=(5.5, 5.5)) 
g = sns.kdeplot(data_to_plot.military_exp, data2=data_to_plot.central_govt_debt) 
g.set_title('Military expenditures & Debt of central governments in 2010')

时间序列数据

前面讨论的横截面数据示例仅来自 2010 年。但是，相反，如果我们仅考虑一个国家（例如美国），并查看其从 2001 年到 2010 年的 10 年间的军费和中央政府债务，那将得到两个时间序列-一个关于美国联邦军费支出以及其他与美国联邦政府债务有关的支出。因此，从本质上讲，时间序列由对单个实体的一个或多个可测量特征的定量观察组成，并在多个时间点获取。在这种情况下，数据代表美国的年度军事支出和政府债务。时间序列数据通常以一些有趣的内部结构为特征，例如趋势，季节性，平稳性，自相关等。这些将在本章的后续部分中进行概念性讨论。

时间序列数据的内部结构需要特殊的表述和分析技术。这些技术将在接下来的章节中通过案例研究和 Python 中工作代码的实现进行介绍。

下图绘制了我们一直在谈论的几个时间序列：

图 1.3：时间序列数据示例

为了生成前面的图，我们将扩展为获取横截面数据图而开发的代码。我们将首先创建两个新的Series，分别代表 1960 年至 2010 年美国军费和中央政府债务的时间序列：

central_govt_debt_us = central_govt_debt.ix[central_govt_debt['Country Code']=='USA', :].T 
military_exp_us = military_exp.ix[military_exp['Country Code']=='USA', :].T

前面代码中创建的两个Series对象合并在一起形成一个DataFrame，并切片以保存 2001 年到 2010 年的数据：

data_us = pd.concat((military_exp_us, central_govt_debt_us), axis=1) 
index0 = np.where(data_us.index=='1960')[0][0] 
index1 = np.where(data_us.index=='2010')[0][0] 
data_us = data_us.iloc[index0:index1+1,:] 
data_us.columns = ['Federal Military Expenditure', 'Debt of Federal  Government'] 
data_us.head(10)

前面的代码准备的数据返回下表：

| | 联邦军事支出 | 联邦政府的债务 |
| 1960 | NaN | NaN |
| 1961 | NaN | NaN |
| 1962 | NaN | NaN |
| 1963 | NaN | NaN |
| 1964 | NaN | NaN |
| 1965 | NaN | NaN |
| 1966 | NaN | NaN |
| 1967 | NaN | NaN |
| 1968 | NaN | NaN |
| 1969 | NaN | NaN |

上表显示，从 1960 年开始的几年中，都无法获得联邦军事开支和联邦债务的数据。因此，在绘制时间序列之前，我们从Dataframe data_us中删除缺少值的行：

data_us.dropna(inplace=True)
print('Shape of data_us:', data_us.shape)

从打印功能的输出中可以看到，在删除缺失值之后，DataFrame 有 23 行：

Shape of data_us: (23, 2)

删除缺少值的行后，我们显示data_us的前十行，如下所示：

| | 联邦军事支出 | 联邦政府的债务 |
| 1988 | 5.57993 | 42.0258 |
| 1989 | 5.37472 | 43.1439 |
| 1990 | 5.12025 | 45.3772 |
| 1991 | 4.53985 | 48.633 |
| 1992 | 4.66626 | 50.6016 |
| 1993 | 4.32693 | 50.1657 |
| 1994 | 3.94129 | 49.3475 |
| 1995 | 3.63849 | 49.2366 |
| 1996 | 3.35074 | 46.7174 |
| 1997 | 3.2099 | 43.2997 |

最后，通过执行以下代码来生成时间序列：

# Two subplots, the axes array is 1-d
f, axarr = plt.subplots(2, sharex=True)
f.set_size_inches(5.5, 5.5)
axarr[0].set_title('Federal Military Expenditure during 1988-2010 (% of GDP)')
data_us['Federal Military Expenditure'].plot(linestyle='-', marker='*', color='b', ax=axarr[0])
axarr[1].set_title('Debt of Federal Government during 1988-2010 (% of GDP)')
data_us['Debt of Federal Government'].plot(linestyle='-', marker='*', color='r', ax=axarr[1])

面板数据

到目前为止，我们已经看到了从多个个人获取的数据，但这些数据是在一个时间点（横截面），或者是从一个个体实体获取的，但数据是在多个时间点（时间序列）。但是，如果我们在多个时间点上观察到多个实体，则会得到面板数据，也称为纵向数据。扩展我们前面有关军事支出的示例，现在让我们考虑 1960-2010 年同期的四个国家。结果数据将是一个面板数据集。下图显示了这种情况下的面板数据。在绘制数据之前，已删除具有缺失值（对应于 1960 年至 1987 年）的行。

图 1.4：面板数据示例

通用面板数据回归模型可以表示为 y_it = W x _it + b + _it，该模型将因变量 y_it 表示为解释变量 x_it，其中 W 是 x_it 的权重，b是偏差项，__it 是误差。i代表针对多个时间点收集数据的个人，这些时间点由j表示。显而易见，这种类型的面板数据分析试图对多个单个时间点和多个时间点的变化进行建模。变化由 _it 反映，假设确定了必要的数学处理方法。例如，如果假设 ϵ_it 相对于i和t非随机变化，则它减小为代表随机噪声的虚拟变量。这种类型的分析称为固定效应模型。另一方面，it_it 随i和t随机变化，需要对误差进行特殊处理，并在随机效应模型中进行处理。

让我们准备绘制上图所需的数据。我们将继续扩展本章前面用于横截面和时间序列数据的代码。我们首先创建一个DataFrame，其中包含上图中提到的四家公司的数据。这样做如下：

chn = data.ix[(data['Indicator Name']=='Military expenditure (% of GDP)')&\
               (data['Country Code']=='CHN'),index0:index1+1
             ]
chn = pd.Series(data=chn.values[0], index=chn.columns)
chn.dropna(inplace=True)

usa = data.ix[(data['Indicator Name']=='Military expenditure (% of GDP)')&\
               (data['Country Code']=='USA'),index0:index1+1
             ]
usa = pd.Series(data=usa.values[0], index=usa.columns)
usa.dropna(inplace=True)

ind = data.ix[(data['Indicator Name']=='Military expenditure (% of GDP)')&\
               (data['Country Code']=='IND'),index0:index1+1
             ]
ind = pd.Series(data=ind.values[0], index=ind.columns)
ind.dropna(inplace=True)

gbr = data.ix[(data['Indicator Name']=='Military expenditure (% of GDP)')&\
               (data['Country Code']=='GBR'),index0:index1+1
             ]
gbr = pd.Series(data=gbr.values[0], index=gbr.columns)
gbr.dropna(inplace=True)

现在已经可以为所有五个国家/地区准备好数据了，我们将使用以下代码对它们进行绘制：

plt.figure(figsize=(5.5, 5.5))
usa.plot(linestyle='-', marker='*', color='b')
chn.plot(linestyle='-', marker='*', color='r')
gbr.plot(linestyle='-', marker='*', color='g')
ind.plot(linestyle='-', marker='*', color='y')
plt.legend(['USA','CHINA','UK','INDIA'], loc=1)
plt.title('Miltitary expenditure of 5 countries over 10 years')
plt.ylabel('Military expenditure (% of GDP)')
plt.xlabel('Years')s

Jupyter 笔记本（具有用于生成所有上述图形的代码）位于 GitHub 存储库中code文件夹下的Chapter_1_Different_Types_of_Data.ipynb。

有关不同类型数据的讨论为进一步研究时间序列奠定了基础。我们将通过了解通常可以在时间序列中找到的数据的特殊属性或其中包含固有时间序列的面板数据来开始这样做。

时间序列的内部结构

在本节中，我们将在概念上解释需要时间序列数据进行特殊数学处理的以下特殊特征：

大势所趋
季节性
周期性运动
意外的变化

大多数时间序列具有一个或多个上述内部结构。基于此概念，时间序列可以表示为 x [t] = f [t] + s [t] + c [t] + e [t]，它是按顺序排列的趋势，季节，周期性和不规则分量的总和。此处，t是在 t = 1,2,3 ... N 连续且等间隔的时间点上进行系列观测的时间索引。

时间序列分析的目的是将时间序列分解为其组成特征，并为每个特征开发数学模型。然后使用这些模型来了解是什么导致了观察到的时间序列行为，并预测了未来时间点的序列。

大势所趋

从长远来看，当一个时间序列显示出向上或向下的移动时，则具有总体趋势。检查总趋势是否存在的快速方法是绘制时间序列，如下图所示，该时间序列显示了从 1974 年到 1987 年测量的空气中 CO [2] 浓度：

图 1.5：CO2 读数的时间序列呈上升趋势

但是，从短期来看，总体趋势可能并不明显。短期影响（例如季节性波动和不规则变化）会导致时间序列重新出现在过去观察到的较低或较高值，因此可能会暂时混淆任何总体趋势。在放大的 1979 年至 1981 年期间，在相同的时间序列中，CO [2] 的浓度很明显，如下图所示。因此，要揭示总体趋势，我们需要一个可以追溯到过去的时间序列。

图 1.6：较短的 CO2 读数时间序列无法显示总体趋势

时间序列的总体趋势是由于它所代表的过程或系统的根本性变化或系统性变化。例如，1974 年至 1987 年期间 CO [2] 浓度的上升趋势可归因于这些年来汽车和工业化的逐步发展。

通常通过将时间序列设置为对时间的回归，并将其他已知因素作为解释变量来对总体趋势进行建模。然后可以将回归线或趋势线用作时间序列的长期运动的预测。趋势线留下的残差将进一步分析其他有趣的属性，例如季节性，周期性行为和不规则变化。

现在，让我们看一下生成关于 CO [2] 浓度的前述图的代码。我们还将展示如何使用线性回归来建立趋势模型，其中以时间指数（在这种情况下为数据中的年份指数）作为解释变量，以 CO [2] 浓度为因变量。但是首先，让我们将数据加载到pandas.DataFrame中。

此示例的数据在 GitHub 存储库的datasets文件夹下的 Excel 文件Monthly_CO2_Concentrations.xlsx中。

我们首先导入所需的软件包，如下所示：

from __future__ import print_function 
import os 
import pandas as pd 
import numpy as np 
%matplotlib inline 
from matplotlib import pyplot as plt 
import seaborn as sns 
os.chdir('D:\Practical Time Series') 
data = pd.read_excel('datasets/Monthly_CO2_Concentrations.xlsx',                                  converters={'Year': np.int32, 'Month':  np.int32}) 
data.head()

我们已将参数converters传递给read_excel函数，以确保为Year和Month列分配了整数（np.int32）数据类型。前面的代码行将生成下表：

| | CO [2] | 年 | 月 |
| 0 | 333.13 | 1974 | 5 |
| 1 | 332.09 | 1974 | 6 |
| 2 | 331.10 | 1974 | 7 |
| 3 | 329.14 | 1974 | 8 |
| 4 | 327.36 | 1974 | 9 |

绘制之前，我们必须删除所有缺少值的列。此外，DataFrame以Year和Month的升序排序。这些操作如下：

data = data.ix[(~pd.isnull(data['CO2']))&\ 
               (~pd.isnull(data['Year']))&\ 
               (~pd.isnull(data['Month']))] 
data.sort_values(['Year', 'Month'], inplace=True)

最后，通过执行以下几行来生成 1974 年至 1987 年时间段的图：

plt.figure(figsize=(5.5, 5.5))
data['CO2'].plot(color='b')
plt.title('Monthly CO2 concentrations')
plt.xlabel('Time')
plt.ylabel('CO2 concentratition')
plt.xticks(rotation=30)

DataFrame之后的这三年是 1980 年至 1981 年期间数据的放大版本：

plt.figure(figsize=(5.5, 5.5))
data['CO2'].loc[(data['Year']==1980) | (data['Year']==1981)].plot(color='b')
plt.title('Monthly CO2 concentrations')
plt.xlabel('Time')
plt.ylabel('CO2 concentratition')
plt.xticks(rotation=30)

接下来，让我们拟合趋势线。为此，我们从scikit-learn导入LinearRegression类，并在时间索引上拟合线性模型：

from sklearn.linear_model import LinearRegression 
trend_model = LinearRegression(normalize=True, fit_intercept=True) 
trend_model.fit(np.array(data.index).reshape((-1,1)), data['CO2']) 
print('Trend model coefficient={} and intercept={}'.format(trend_model.coef_[0], 
                                                           trend_model.intercept_) 
      )

这将产生以下输出：

Trend model coefficient=0.111822078545 and intercept=329.455422234

下图显示了从趋势线模型获得的残差，这些残差似乎具有季节性行为，这将在下一部分讨论。

残差是通过下面的代码行在前面计算和绘制的：

residuals = np.array(data['CO2']) - trend_model.predict(np.array(data.index).reshape((-1,1)))
plt.figure(figsize=(5.5, 5.5))
pd.Series(data=residuals, index=data.index).plot(color='b')
plt.title('Residuals of trend model for CO2 concentrations')
plt.xlabel('Time')
plt.ylabel('CO2 concentratition')
plt.xticks(rotation=30)

图 1.7：来自二氧化碳读数总体趋势的线性模型中的残差

季节性

季节性表现为时间序列中的重复和周期变化。在大多数情况下，探索性数据分析揭示了季节性的存在。让我们回顾一下 CO [2] 浓度的趋势下降的时间序列。尽管去趋势线系列具有恒定的均值和恒定的方差，但是它以可预测的方式系统地偏离趋势模型。

季节性表现为周期性偏差，例如在 CO [2] 排放量下降趋势观察中看到的那些。圣诞礼物或季节性服装等季节性商品的月销售量的高峰和低谷是具有季节性的时间序列的另一个示例。

确定季节性的一种实用技术是通过以下图进行探索性数据分析：

运行序列图
季节性子系列剧情
多箱图

运行序列图

一个简单的原始时间序列的运行序列图，其中时间在x轴上，变量在y轴上，可以很好地指示时间序列的以下属性：

系列平均运动
方差变化
异常值的存在

下图是从数学公式 x [t] =（At + B）sin（t）+Є（t）获得的假设时间序列的运行序列图时间相关的均值和误差 Є（t），它随正态分布 N（0，at + b）的变化而变化。此外，还包括一些异常高和低的观察值作为异常值。

在这种情况下，运行序列图是识别序列均值和方差以及离群值的有效方法。 CO2 浓度的下降趋势时间序列图是运行序列图的一个示例。

季节性子系列剧情

对于已知的季节性变化周期性，季节性子系列会在连续时间段的批次中重绘原始系列。例如，CO [2] 浓度的周期为 12 个月，因此，下图显示了残差均值和标准差的季节性子系列图。为了可视化残差的季节性，我们创建季度均值和标准差。

一个季节性子系列揭示了两个属性：

连续几个月内，季节内的变化
连续几批之间的季节之间的变化

图 1.8：二氧化碳读数总体趋势线性模型中残差的季度均值

图 1.9：残留物的季度标准偏差与二氧化碳读数总体趋势的线性模型

现在让我们描述用于生成前面的图的代码。首先，我们需要将残留物和四分之一标记添加到 CO2 浓度DataFrame中。这样做如下：

data['Residuals'] = residuals 
month_quarter_map = {1: 'Q1', 2: 'Q1', 3: 'Q1', 
                     4: 'Q2', 5: 'Q2', 6: 'Q2', 
                     7: 'Q3', 8: 'Q3', 9: 'Q3', 
                     10: 'Q4', 11: 'Q4', 12: 'Q4' 
                    } 
data['Quarter'] = data['Month'].map(lambda m: month_quarter_map.get(m))

接下来，通过对Year和Quarter上的数据进行分组来计算季节性平均值和标准差：

seasonal_sub_series_data = data.groupby(by=['Year', 'Quarter'])['Residuals']\ 
                           .aggregate([np.mean, np.std])

这将创建新的DataFrame作为seasonal_sub_series_data，这些年份具有季度平均值和标准偏差。这些列的重命名如下：

seasonal_sub_series_data.columns = ['Quarterly Mean', 'Quarterly Standard Deviation']

接下来，通过运行以下代码行来绘制季度均值和标准差：

#plot quarterly mean of residuals
plt.figure(figsize=(5.5, 5.5))
seasonal_sub_series_data['Quarterly Mean'].plot(color='b')
plt.title('Quarterly Mean of Residuals')
plt.xlabel('Time')
plt.ylabel('CO2 concentratition')
plt.xticks(rotation=30)

#plot quarterly standard deviation of residuals
plt.figure(figsize=(5.5, 5.5))
seasonal_sub_series_data['Quarterly Standard Deviation'].plot(color='b')
plt.title('Quarterly Standard Deviation of Residuals')
plt.xlabel('Time')
plt.ylabel('CO2 concentratition')
plt.xticks(rotation=30)

多箱图

使用季节性箱形图重绘时，季节性子系列图可以提供更多信息，如下图所示。箱形图显示了一批时间单位内季节性数据内的集中趋势和离散度。此外，两个相邻箱形图之间的分离揭示了季节内的变化：

图 1.10：来自二氧化碳读数的一般趋势的线性模型的差值差值

生成箱形图的代码如下：

plt.figure(figsize=(5.5, 5.5)) 
g = sns.boxplot(data=data, y='Residuals', x='Quarter') 
g.set_title('Quarterly Mean of Residuals') 
g.set_xlabel('Time') 
g.set_ylabel('CO2 concentratition')

周期性变化

周期性变化是每隔几个时间单位观察到的运动，但其发生频率低于季节性波动。与季节性不同，周期性变化可能没有固定的变化周期。此外，周期性变化的平均周期会更大（最常见的年份是几年），而同一年内观察到季节性变化，并且对应于时间的年度划分，例如季节，季度以及节日和节假日等。

需要一个时间序列的长期图来识别可能发生的周期性变化，例如，每隔几年出现一次，并表现为重复的波峰和波谷。在这方面，与经济和商业有关的时间序列通常会显示出与通常的商业和宏观经济周期相对应的周期性变化，例如经济衰退期，随后的繁荣时期，但相隔数年的时间跨度。与一般趋势类似，识别周期性运动可能需要追溯到过去的数据。

下图说明了 1960 年至 2016 年期间印度和美国的消费者价格指数（CPI）通货膨胀发生的周期性变化。这两个国家的 CPI 通货膨胀率均表现出周期性变化，大约需要 2-2.5 年的时间。此外，印度的 CPI 通货膨胀率在 1990 年前比 1990 年以后更大。

图 1.11：时间序列数据中的周期性运动示例

来源：上图的数据已从这个页面下载，该网站维护了来自多个主题的时间序列数据。

您可以在 GitHub 存储库中datasets文件夹中的inflation-consumer-prices-annual.xlsx文件中找到 CPI 通胀数据集。

生成该图所编写的代码如下：

inflation = pd.read_excel('datasets/inflation-consumer-prices-annual.xlsx', parse_dates=['Year']) 
plt.figure(figsize=(5.5, 5.5)) 
plt.plot(range(1960,2017), inflation['India'], linestyle='-', marker='*', color='r') 
plt.plot(range(1960,2017), inflation['United States'], linestyle='-', marker='.', color='b') 
plt.legend(['India','United States'], loc=2) 
plt.title('Inflation in Consumer Price Index') 
plt.ylabel('Inflation') 
plt.xlabel('Years')

意外的变化

关于将时间序列表示为四个成分之和的模型，值得注意的是，尽管能够考虑其他三个成分，但我们可能仍会留下不可约的误差成分，该成分是随机的并且不会表现出对时间指标的系统依赖性。这第四部分反映了时间序列中的意外变化。出乎意料的变化是随机的，无法在数学模型中进行明确的未来预测。这种类型的错误是由于缺乏有关可以对这些变化进行建模的解释变量的信息，或者是由于存在随机噪声。

本节中开发的所有代码的 IPython 笔记本都在本书 GitHub 存储库的code文件夹中的Chapter_1_Internal_Structures.ipynb文件中。

时间序列分析模型

时间序列分析的目的是开发一个数学模型，该模型可以解释观察到的时间序列行为，并可能预测该序列的未来状态。所选模型应能够说明可能存在的一个或多个内部结构。为此，我们将概述以下通常用作时间序列分析构建块的常规模型：

零均值模型
随机漫步
趋势模型
季节性模型

零均值模型

零均值模型具有恒定的均值和恒定的方差，并且没有可预测的趋势或季节性。假设来自零均值模型的观测值是独立且均匀分布的（iid），并且表示固定均值周围的随机噪声，该噪声已从时间序列中扣除为常数学期。

让我们考虑 X [1] ，X [2] ，...，X [n] 代表与 a 的 n 个观测值相对应的随机变量。零均值模型。如果 x [1] ，x [2] ，...，x [n] 是来自零平均时间序列的 n 个观测值，则联合观测值的分布是每个时间索引的概率质量函数的乘积，如下所示：

**P(X1 = x1,X2 = x2 , ... , Xn = xn) = f(X1 = x1) f(X2 = x2) ... f(Xn = xn)**

最常见的 f（X [t] = x [t] ）由均值零和方差的正态分布建模σ ²，被认为是模型的不可减少的误差，因此被视为随机噪声。下图显示了零均值序列的方差的正态分布随机噪声：

图 1.12：零均值时间序列

上图是通过以下代码生成的：

import os 
import numpy as np 
%matplotlib inline 
from matplotlib import pyplot as plt 
import seaborn as sns 
os.chdir('D:/Practical Time Series/') 
zero_mean_series = np.random.normal(loc=0.0, scale=1., size=100)

具有恒定方差的零均值表示随机噪声，该随机噪声可以假定无限可能的实数值，并且适合于表示连续变量的时间序列中的不规则变化。但是，在许多情况下，系统或过程的可观察状态本质上可能是离散的，并且被限制为有限数量的可能值 s [1] s [2] 等。。，s [m]。在这种情况下，假设观测变量（X）服从多项式分布，P（X = s [1] ）= p [1] ， P（X = s [2] ）= p [2] ，…，P（X = s [m] ）= p [m] 这样 p [1] + p [2] + ... + p [m] = 1。这样的时间序列是离散的随机过程。

随时间推移多次掷骰子是离散随机过程的一个示例，任何掷骰都有六个可能的结果。一个更简单的离散随机过程是一个二进制过程，例如扔硬币，例如只有两个结果，即正面和反面。下图显示了投掷有偏骰子的模拟过程的 100 次运行，为此，举起偶数面的几率要高于显示奇数面的几率。请注意，平均而言，与奇数面孔的出现次数相比，偶数面孔的发生次数更高。

随机漫步

随机游走作为 n 个 id 的总和给出，它具有零均值和恒定方差。基于此定义，时间总和 t 处的随机游走的实现由总和 S = x [1] + x [2] + ... + x [n]。下图显示了从 iid 获得的随机游动，其随零均值和单位方差的正态分布而变化。

随机游走很重要，因为如果在一个时间序列中发现这种行为，则可以通过将来自两个连续时间指标的观测值的差异作为 S [t] -S 轻松地简化为零均值模型 [t-1] = x [t] 是具有零均值和恒定方差的 iid。

图 1.13：随机游走时间序列

上图中的随机游走可以通过采用上一节中讨论的零均值模型的累加和来生成。以下代码实现了这一点：

random_walk = np.cumsum(zero_mean_series) 
plt.figure(figsize=(5.5, 5.5)) 
g = sns.tsplot(random_walk) 
g.set_title('Random Walk') 
g.set_xlabel('Time index')

趋势模型

这种类型的模型旨在捕获时间序列中的长期趋势，可以将其拟合为时间指数的线性回归。当时间序列没有任何周期性或季节性波动时，可以将其表示为趋势和零均值模型的总和，表示为 x [t] =μ（t）+ y [t]，其中μ（t）是该序列的时间依赖性长期趋势。

趋势模型μ（t）的选择对于正确捕获时间序列的行为至关重要。探索性数据分析通常为假设模型在t中是线性还是非线性提供了提示。线性模型为μ（t）= wt + b，而二次模型为μ（t）= w [1] t + w [2] t ² + b。有时，可以根据时间指数，例如μ（t）= w [0] t ^p + b，通过更复杂的关系来推测趋势。

通过以t为解释变量，以μ作为解释变量进行回归，可以获得趋势模式（如先前讨论的趋势）中的权重和偏差。趋势模型的残差 x [t] -μ（t）被视为不可归约的噪声，并且作为零均值分量 y [ [t]。

季节性模型

季节性表现为时间序列中的周期性和重复性波动，因此可以将其建模为已知周期性的正弦波的加权和。假设长期趋势已被趋势线消除，则季节性模型可以表示为x[t] = s[t] + y[t]，其中具有已知周期性的季节性变化为α：

季节性模型也称为谐波回归模型，因为它们试图拟合多个正弦波的总和。

此处描述的四个模型是成熟的时间序列模型的构建基块。如您现在可能已经收集到的，零和模型代表了系统的不可减少的误差，所有其他三个模型旨在通过适当的数学转换将给定的时间序列转换为零和模型。为了获得有关原始时间序列的预测，应用了相关的逆变换。

接下来的章节将详细介绍此处讨论的四种模型。但是，我们可以在以下四个步骤中总结时间序列分析的通用方法：

可视化时间索引不同粒度的数据，以揭示长期趋势和季节性波动
拟合趋势线可捕获长期趋势并绘制残差以检查季节性或不可减少的误差
拟合谐波回归模型以捕获季节性
绘制季节性模型剩余的残差以检查不可减少的误差

这些步骤通常足以为大多数时间序列开发数学模型。各个趋势和季节性模型可以是简单的也可以是复杂的，具体取决于原始时间序列和应用。

可在本书 GitHub 存储库的code文件夹中的Chapter_1_Models_for_Time_Series_Analysis.ipynb IPython 笔记本中找到本节中编写的代码。

自相关和部分自相关

在应用了上一节中讨论的数学变换之后，我们经常会得到平稳（或弱平稳）时间序列，其特征是恒定的均值 E（x [t] ）和相关性仅取决于两个时间步长之间的时滞，而与时间步长的值无关。这种协方差是时间序列分析的关键，当归一化为-1 到 1 时，称为自协方差或自相关，因此，自相关表示为二阶矩 E（x [t] ，x [t + h] ）= g（h）显然仅是时间滞后h的函数，并且与实际时间指标t无关。自相关的这种特殊定义确保了它是与时间无关的属性，因此可以可靠地用于推断时间序列的未来实现。

自相关反映了索引为t的时间序列与索引为 t-h 或 t + h 的时间序列之间的线性相关程度。正自相关表示时间序列的当前值和将来值沿相同方向移动，而负值表示时间序列的当前值和将来值沿相反方向移动。如果自相关接近于零，则可能很难找到序列中的时间相关性。由于此属性，自相关在预测h时间提前的时间序列的未来状态时很有用。

通过绘制给定时间序列的自相关函数（ACF）的观察值，可以确定自相关的存在。该图通常称为 ACF 图。让我们说明绘制 ACF 的观察值如何帮助检测自相关的存在。为此，我们首先绘制 2016 年 1 月至 2016 年 12 月观察到的道琼斯工业平均指数（DJIA）的每日价值：

图 1.14：道琼斯工业平均指数的时间序列

从上图可以明显看出，道琼斯工业平均指数开始上升时，会持续上升一段时间，反之亦然。但是，我们必须通过 ACF 图来确定这一点。

该图的数据集已从 Yahoo Finance 下载，并以DJIA_Jan2016_Dec2016.xlsx的形式保存在本书 GitHub 存储库的datasets文件夹下。

我们将使用pandas从 Excel 文件中读取数据，并与 matplotlib 一起使用 seaborn 来可视化时间序列。像以前一样，我们还将使用 os 包来设置工作目录。

因此，让我们首先导入以下软件包：

import os 
import pandas as pd 
%matplotlib inline 
from matplotlib import pyplot as plt 
import seaborn as sns  
os.chdir('D:/Practical Time Series')

接下来，我们将数据加载为pandas.DataFrame并显示其前 10 行以查看数据集的列：

djia_df = pd.read_excel('datasets/DJIA_Jan2016_Dec2016.xlsx') 
djia_df.head(10)

上面的代码显示数据集的前 10 行，如下表所示：

| | 日期 | 打开 | 高 | 低 | 关闭 | 调整关闭 | 音量 |
| 0 | 2016-01-04 | 17405.480469 | 17405.480469 | 16957.630859 | 17148.939453 | 17148.939453 | 148060000 |
| 1 | 2016-01-05 | 17147.500000 | 17195.839844 | 17038.609375 | 17158.660156 | 17158.660156 | 105750000 |
| 2 | 2016-01-06 | 17154.830078 | 17154.830078 | 16817.619141 | 16906.509766 | 16906.509766 | 120250000 |
| 3 | 2016-01-07 | 16888.359375 | 16888.359375 | 16463.630859 | 16514.099609 | 16514.099609 | 176240000 |
| 4 | 2016-01-08 | 16519.169922 | 16651.890625 | 16314.570313 | 16346.450195 | 16346.450195 | 141850000 |
| 5 | 2016-01-11 | 16358.709961 | 16461.849609 | 16232.030273 | 16398.570313 | 16398.570313 | 127790000 |
| 6 | 2016-01-12 | 16419.109375 | 16591.349609 | 16322.070313 | 16516.220703 | 16516.220703 | 117480000 |
| 7 | 2016-01-13 | 16526.630859 | 16593.509766 | 16123.200195 | 16151.410156 | 16151.410156 | 153530000 |
| 8 | 2016-01-14 | 16159.009766 | 16482.050781 | 16075.120117 | 16379.049805 | 16379.049805 | 158830000 |
| 9 | 2016-01-15 | 16354.330078 | 16354.330078 | 15842.110352 | 15988.080078 | 15988.080078 | 239210000 |

上表中的第一列始终是pandas.read_csv函数创建的默认行索引。

我们使用Close列中给出的 DJIA 的关闭值来说明自相关和 ACF 函数。时间序列图已生成如下：

plt.figure(figsize=(5.5, 5.5)) 
g = sns.tsplot(djia_df['Close']) 
g.set_title('Dow Jones Industrial Average between Jan 2016 - Dec 2016') 
g.set_xlabel('Time') 
g.set_ylabel('Closing Value')

接下来，通过为滞后h的不同值计算自相关来估计 ACF，在这种情况下，滞后h在 0 到 30 之间变化。Pandas.Series.autocorr函数用于为滞后不同值计算自相关。给出的代码如下：

lag = range(0,31) 
   djia_acf = [] 
for l in lag: 
    djia_acf.append(djia_df['Close'].autocorr(l))

前面的代码在从 0 到 30 的 31 个滞后值列表中进行迭代。滞后 0 表示观察值与自身之间具有自相关性（换言之为自相关），因此也可以确定为 1.0 在下图中。 DJIA 收盘价的自相关似乎随着延迟而线性下降，并且下降速度在 18 天左右出现明显变化。在 30 天的滞后时间内，ACF 略高于 0.65。

图 1.15：针对各种滞后计算出的道琼斯工业平均指数的自相关

前面的图是通过以下代码生成的：

plt.figure(figsize=(5.5, 5.5)) 
g = sns.pointplot(x=lag, y=djia_acf, markers='.') 
g.set_title('Autocorrelation function for DJIA') 
g.set_xlabel('Lag in terms of number of trading days') 
g.set_ylabel('Autocorrelation function') 
g.set_xticklabels(lag, rotation=90) 
plt.savefig('plots/ch2/B07887_02_11.png', format='png', dpi=300)

ACF 图显示，在 DJIA Close 值的情况下，自相关与观测之间的时间间隔具有函数相关性。

为运行本节中的分析而开发的代码位于本书 GitHub 存储库的code文件夹下的 IPyton 笔记本Chapter1_Autocorrelation.ipynb中。

我们编写了一个 for 循环来计算不同滞后的自相关，并使用sns.pointplot函数绘制结果。或者，statsmodels.graphics.tsaplots的plot_acf函数可计算和绘制各种滞后的自相关。此外，该功能还可以绘制 95％的置信区间。在这些置信区间之外的自相关在统计上是显着的相关，而在置信区间内的自相关是由于随机噪声引起的。
由plot_acf生成的自相关和置信区间如下图所示：

图 1.16：95％置信区间的道琼斯工业平均指数自相关

到目前为止，我们已经讨论了自相关，它是变量 x_t 和 x_（t + h）之间线性相关性的量度。 自回归（AR）模型以 x_（t + h）和 x_t 之间的线性回归的形式捕获了这种依赖性。但是，时间序列倾向于逐步携带信息和相关性结构，因此滞后h的自相关也受到中间变量 x_t，x_（t + 1）…x_（t + h- 1）。因此，在存在中间变量的情况下，自相关不是 x_t 和 x_（t + h）之间相互关系的正确度量。因此，基于自相关在 AR 模型中选择h是错误的。当消除了中间变量的影响时，部分自相关通过测量 x_t 和 x_（t + h）之间的相关性来解决此问题。因此，时间序列分析中的部分自相关定义了 x_t 和 x_（t + h）之间的相关性，但并未由滞后 t + 1 到 ] t + h-1。
局部自相关有助于确定 AR（h）模型的阶数 h。让我们使用 plot_pacf 绘制 DJIA 关闭值的局部自相关，如下所示：

图 1.17：道琼斯工业平均指数与 95％置信区间的部分自相关

滞后零处的第一个部分自相关始终为 1.0。从上图中可以看出，仅在滞后一个时，部分自相关在统计上就很显着，而对于其他滞后，则在 95％的置信区间内。因此，对于 DJIA Close Values，AR 模型的顺序为 1。

概括

在本章中，我们讨论了几种数据类型，例如横截面，时间序列和面板数据。我们深入研究了使时间序列数据变得特别的特殊属性。已经讨论了 Python 中的几个示例和工作代码，以了解如何在时间序列上进行探索性数据分析以可视化其属性。我们还描述了自相关和部分自相关以及图形技术来检测时间序列中的这些。本章讨论的主题为我们提供了一个使用 Python 进行时间序列数据更详细讨论的阶段。
在下一章中，您将学习如何按时间序列读取更复杂的数据类型，以及如何使用此类信息进行更深入的探索性数据分析。我们将在时间序列平稳的情况下重新审视自相关。将讨论检测自相关的统计方法。我们还将讨论平稳性的重要性，并描述用于平稳化非平稳时间序列的不同微分和平均方法。讨论了时间分解的加法和乘法模型，用于估计趋势和季节性。

二、了解时间序列数据

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

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

生成固定数据对于增强时间序列预测模型非常重要。趋势，季节和周期性成分的推论将使我们面临不规则的波动，而不能仅以时间指数作为解释变量来建模。因此，为了进一步改善预测，不规则波动被假定为独立于且分布均匀的（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：减去季度啤酒产量与其移动平均值之差后剩余残留物的自相关

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

Take seasonal 2 × F_hat[t]^(4) MA.
通过从原始时间序列中删除2 × F_hat[t]^(4) MA 来计算残差。
检查残差的 ACF 的随机性。
如果残差的 ACF 已经是随机的，则残差是固定的，如果没有，则继续下一步。
取周期性为 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 进行加性和乘性模型的时间序列分解。
在下一章中，我们将深入研究处理噪声时间序列数据的指数平滑技术。

三、基于指数平滑的方法

本章介绍了对时间序列信号的数据平滑处理。本章的结构如下：

时间序列平滑介绍
一阶指数平滑
二阶指数平滑
建模高阶指数平滑
概括

时间序列平滑介绍

时间序列数据由信号和噪声组成，其中信号捕获了过程的内在动力。但是，噪声代表信号的未建模分量。时间序列信号的内在动力学可以与过程的平均值一样简单，也可以是观测值内的复杂功能形式，如下所示：

*x[t] = f(x[i]) + ε[t]* for *i=1,2,3, ... t-1*

此处，x [t] 是观察值，而ε [t] 是白噪声。f（x [i] ）表示功能形式；常量作为函数形式的示例如下：

x [t] =μ+ε [t]

此处，上式中的常数值μ用作漂移参数，如下图所示：

图 3.1：带有漂移参数的时间序列示例

由于ε [t] 是白噪声，因此这种基于平滑的方法可以通过消除固有函数形式与随机噪声进行分离。平滑预测方法可以视为采用输入并分离趋势和噪声成分的过滤器，如下图所示：

图 3.2：平滑过程的框架

估计趋势和噪声的提取效率取决于与时间序列信号组成相关的其他参数，例如趋势，季节性和残差（噪声）的存在。为了处理时间序列的每个组成部分，需要进行不同的处理。本章将介绍多种平滑处理时间序列信号不同成分的方法。使用从 1946 年 1 月到 1959 年 12 月的每月水平收集的纽约出生数据集，演示了一个由趋势，季节性和白噪声组成的时间序列信号示例：

Import requests 
import statsmodels.api as sm 
import io 
import pandas as pd 

# Load Dataset 
DATA_URL="http://robjhyndman.com/tsdldata/data/nybirths.dat" 
fopen = requests.get(DATA_URL).content 
ds=pd.read_csv(io.StringIO(fopen.decode('utf-8')),  header=None, names=['birthcount']) 
print(ds.head()) 

# Add time index 
date=pd.date_range("1946-01-01", "1959-12-31", freq="1M") 
ds['Date']=pd.DataFrame(date) 
ds = ds.set_index('Date') 

# decompose dataset 
res = sm.tsa.seasonal_decompose(ds.birthcount, model="additive") 
resplot = res.plot()

代码中的不同功能如下使用：

前面脚本中的requests.get函数用于从 URL DATA_URL 获取数据。
要处理数据集，请使用 pandas DataFrame。
统计模型模块的seasonal_decompose功能用于将时间序列信号分解为趋势，季节性和残差成分。如第 2 章，“时间序列简介”中所述，分解可以是加法运算或乘法运算。下图显示了信号不同成分的示例：

图 3.3：时间序列信号的分解

先前的时间序列信号由趋势，季节性和残差（噪声）组成。如图进行平滑处理的框架所示，平滑有助于去除残留分量，并捕获预测信号的趋势和季节性分量。

包括均值，趋势和非季节性模式的第一步模型是使用平滑外推法。第 2 章，“了解时间序列数据”中讨论了使用移动平均的基本平滑。移动平均平滑使用先前的所有观察评估期望 E（xt），如下所示：

通常，在窗口上执行简单的移动平均；因此，使用目标函数在最佳窗口上评估估计的预测，以最大程度地减少误差：

平滑方法基于以下假设：时间序列数据是局部平稳的，且均值变化很小。因此，我们可以使用时间t的平均值来预测时间[delta]小到足以保持信号平稳的 t + 1。该模型是没有漂移模型的均值和随机游走之间的折衷方案。该模型也称为平滑模型，因为它可以平滑数据集中的颠簸。这些基于移动平均的方法的主要局限性在于，它会同等地对待平滑中使用的所有n样本，而忽略了观察新近度的影响，即，将较高的权重赋给最近的观察，如以下等式所示：

w [1] > w [2] > w [3] > ... > w [T [] 和T是窗口大小。这对权重的评估提出了另一个挑战。通过对观察值使用指数衰减权重进行指数平滑，可以解决移动平均值和加权移动平均值的局限性。

一阶指数平滑

一阶指数平滑或简单指数平滑适合具有恒定方差且无季节性的情况。通常建议使用此方法进行短期预测。第 2 章和“了解时间序列数据”引入了朴素的预测方法，其中将地平线 h 的预测定义为 t 值（或最后一个观测值）：

*x[t+h] = x[t]*

该方法通过简单的移动平均值进行扩展，该方法使用多个历史点的均值扩展了朴素的方法：

该方法假定所有历史观测值都具有相同的权重，如下图所示：

图 3.4：随窗口大小增加分配给观察的权重

随着移动平均的窗口大小的增加，分配给每个观察值的权重会变小。一阶指数通过为历史数据点提供指数，从而扩展了当前方法，历史数据点的权重从最新数据点到最早的数据点呈指数递减。一阶指数平滑可以定义如下：

α是[0,1]之间的平滑因子，它控制权重降低的速率，x [t] 是时间时的观测值 ] t。下图显示了平滑过程，它显示了具有不同平滑因子α的权重的衰减：

图 3.5：以不同的 alpha 值分配给历史观测值的权重示意图

α的值越高，权重衰减越快；因此，历史数据对预测值的影响较小。单阶指数平滑的预测方程式可以进一步简化为：

*F[1] = x[1]*

*F[2] = αx[1] + (1-α)F[1] = αx[1] + (1-α)x[1] = x[1]*

*F[3] = αx[2] + (1-α)F[2]*

...

*F[t-1] = αx[t] + (1-α)F[t-1]*

在此，F [t-1] 是在时间t的预测值，表示如下：

简化形式显示了时间 t 的预测值作为时间 t-1 处的函数预测值与时间 t 处的观测值之间的关系。该预测也称为 Holt-Winters 预测。

从前面的等式中可以看出，在第一时间实例中，将初始值分配给第一预测值作为观察值，F [1] = x [1]，第二个时间实例也采用第一个实例的值。重新引导预测的另一种方法是使用可用数据或数据子集的平均值，如下所示：

在移动平均指数平滑中检测窗口大小和权重的类似问题要求优化平滑参数α。正确选择α非常关键。例如，选择α= 0 会将所有值平滑为初始分配的值，如下所示：

*F[t+1] = F[t-1]*

类似地，α= 1 表示指数平滑的最平滑情况，如下所示：

*F[t+1] = x[t]*

让我们以一个 IBM 普通股收盘示例为例，使用单一平滑方法进行预测。第一步是加载所需的模块：

# Load modules 
from __future__ import print_function 
import os 
import pandas as pd 
import numpy as np 
from matplotlib import pyplot as plt

当前示例将使用 Python 中的os，pandas，numpy和matplotlib模块。使用pandasDataFrame将数据集加载到 Python 环境中：

# Load Dataset 
ibm_df = pd.read_csv('datasets/ibm-common-stock-closing-prices.csv') 
ibm_df.head()

为了方便起见，我们将重命名这些列：

#Rename the second column 
ibm_df.rename(columns={'IBM common stock closing prices': 'Close_Price'},inplace=True)

输出数据集将如图图 3.6 所示：

图 3.6：样本 IBM 收盘股价数据集

在单平滑指数方法中，预测值生成如下：

...

前面的系列可以通过以下方式在 Python 中实现：

# Function for Single exponential smoothing 
def single_exp_smoothing(x, alpha): 
    F = [x[0]] # first value is same as series 
    for t in range(1, len(x)): 
        F.append(alpha * x[t] + (1 - alpha) * F[t-1]) 
    return F

设置有初始预测值的single_exp_smoothing功能被指定为该系列的第一个值。让我们首先评估α= 0 和α= 1 的极端预测情况：

# Single exponential smoothing forecasting 
ibm_df['SES0'] = single_exp_smoothing(ibm_df['Close_Price'], 0) 
ibm_df['SES0'] = single_exp_smoothing(ibm_df['Close_Price'], 1)

上一个脚本的输出如下：

图 3.7：平滑参数为零和一的简单指数平滑

上图说明，在α= 0 时，预测值是一个常数，而对于α= 1，预测的序列偏移1时滞。。可以对平滑值 0.2 的单个平滑预测进行如下评估：

# Single exponential smoothing forecasting 
ibm_df['SES'] = single_exp_smoothing(ibm_df['Close_Price'], 0.8)

可以针对实际数据绘制指数平滑的结果，如图图 3.18 所示：

### Plot Single Exponential Smoothing forecasted value 
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['SES'].plot(ax=ax, color='r') 
ax.set_title('Single Exponential Smoothing')

图 3.8：实际值与预测值之间的比较

可以使用标准目标函数，例如平均平方误差（MSE 或平均绝对误差（MAS）：

同样，MAS 的评估如下：

让我们评估α对拟合的影响。为了评估它，使用不同的平滑参数开发了多个模型，如下所示：

#Calculate the single exponential forecast at different values 
ibm_df['SES2']  = single_exp_smoothing(ibm_df['Close_Price'], 0.2) 
ibm_df['SES6']= single_exp_smoothing(ibm_df['Close_Price'], 0.6) 
ibm_df['SES8']= single_exp_smoothing(ibm_df['Close_Price'], 0.8)

图 3.5 将预测值与实际值进行比较：

图 3.9：alpha 的效果示意图

上图说明了 alpha 对预测的影响很大；因此，在设置预测时获取正确的 alpha 值至关重要。可以使用以下脚本获得以上图表：

# Plot the curves 
f, axarr = plt.subplots(3, sharex=True) 
f.set_size_inches(5.5, 5.5) 

ibm_df['Close_Price'].iloc[:45].plot(color='b', linestyle = '-', ax=axarr[0]) 
ibm_df['SES2'].iloc[:45].plot(color='r', linestyle = '-', ax=axarr[0]) 
axarr[0].set_title('Alpha 0.2') 

ibm_df['Close_Price'].iloc[:45].plot(color='b', linestyle = '-', ax=axarr[1]) 
ibm_df['SES6'].iloc[:45].plot(color='r', linestyle = '-', ax=axarr[1]) 
axarr[1].set_title('Alpha 0.6') 

ibm_df['Close_Price'].iloc[:45].plot(color='b', linestyle = '-', ax=axarr[2]) 
ibm_df['SES8'].iloc[:45].plot(color='r', linestyle = '-', ax=axarr[2]) 
axarr[2].set_title('Alpha 0.8')

由于平滑有助于减少数据集方差，因此它将减少预测序列在零到实际数据集方差之间的方差：

解决前面的方程式将导致以下差异：

在此，T是时间序列的长度。对于系列 x [T] 的单位方差，由预测系列捕获的方差将基于平滑参数α改变，如下所示：

图 3.10：使用 alpha 变化的简单指数平滑捕获的方差

二阶指数平滑

如果一阶指数平滑效果不佳，则时间序列数据中存在趋势。这种趋势通常在许多领域都可以观察到，例如当电子商务公司进行营销活动时，销售的增长或公司的任何良好的年度业绩都会对其股票价格产生看涨的影响。线性趋势可能由于时间和响应之间的线性趋势而发生：

*x[t] = constant + ωt + ε[t]*

在此，ω是导致趋势的系数。通过将另一项包括到一阶指数平滑中，二阶指数平滑可以帮助捕获时间序列数据的趋势，如下所示：

在此，T [t] 捕获指数平滑的趋势分量，并表示如下：

在此，α是数据平滑因子，β是趋势平滑因子，其值在[0，1]之间。下一阶段的预测可以如下生成：

在二阶平滑中，趋势分量的初始值可以通过多种方式分配：

在此，n是观察次数。让我们演示一个二阶平滑的示例。我们将使用啤酒产量数据。以下是演示步骤：

第一步是加载所需的模块：

    # Load modules 
    from __future__ import print_function 
    import os 
    import pandas as pd 
    import numpy as np 
    from matplotlib import pyplot as plt 

    Load beer dataset as pandas DataFrame. 
    #Read dataset into a pandas.DataFrame 
    beer_df = pd.read_csv('datasets/quarterly-beer-production-in-
    aus-March 1956-June 1994.csv')

让我们创建一个用于双指数平滑的函数：

    # Function for double exponential smoothing 
    def double_exp_smoothing(x, alpha, beta): 
        yhat = [x[0]] # first value is same as series 
        for t in range(1, len(x)): 
            if t==1: 
                F, T= x[0], x[1] - x[0] 
            F_n_1, F = F, alpha*x[t] + (1-alpha)*(F+T) 
            T=beta*(F-F_n_1)+(1-beta)*T 
            yhat.append(F+T) 
        return yhat

前面的函数将时间序列作为带有 alpha 和 beta 的输入。前面的实现使用前两次出现的差来设置初始趋势值。

让我们评估一下 alpha beta 边界情况的性能，即 alpha 和 beta 平滑参数的（0,0），（0，1），（1,0）和（1,1）值：

    # Effect of alpha and beta 
    beer_df['DEF00'] = double_exp_smoothing(beer_df['Beer_Prod'],0, 0) 
    beer_df['DEF01'] = double_exp_smoothing(beer_df['Beer_Prod'],0, 1) 
    beer_df['DEF10'] = double_exp_smoothing(beer_df['Beer_Prod'],1, 0) 
    beer_df['DEF11'] = double_exp_smoothing(beer_df['Beer_Prod'],1, 1)

下图显示了适合二阶指数平滑的前一个脚本的结果：

图 3.11：二阶指数平滑中的 alpha 和 beta 平滑参数对啤酒销售数据集的影响

当α= 0 时，初始值保持恒定；因此，趋势参数不起作用。但是，当α= 1 和β= 1 时，二阶指数平滑可以写为：

时间t时的预测取决于先前的值和趋势成分。当β设置为零时，t-1 的趋势分量将取决于 t-2：

因此，趋势分量值取决于初始分配值，并且是一个常数。

类似地，对于α= 1 和β= 1，二阶指数平滑可以简化如下：

与α的配置相比，在α= 1 和β= 1 的情况下，将 t-2 和 t-3 的差与时间t预测值相加。 α= 1 和β= 0，这使预测值发生偏移，并且更接近于实际预测值。

让我们使用α和β的中间值对啤酒数据执行双指数平滑，如下所示：

beer_df['DEF'] = double_exp_smoothing(beer_df['Beer_Prod'], 0.4, 0.7)

图 3.12：实际和双指数预测值之间的比较

我们还比较一下单指数和双指数平滑的性能：

图 3.13：单指数和双指数预测值之间的比较

上图显示，与单指数平滑相比，双指数平滑能够更好地捕获当前数据集的真实信号变化。但是，在趋势成分趋于零的情况下，单指数和双指数平滑方法的性能是可比的。

建模高阶指数平滑

该概念可以进一步扩展为使用n^th 阶多项式模型的高阶指数平滑：

在此，误差ε [t] 〜N（0，σ ² ）正态分布，均值为 0，σ ^[2] 方差。用于更高阶的指数平滑器如下：

...

这是平滑器的权重。通常，不及时使用高阶指数平滑，因为即使对于二阶平滑，计算也变得非常困难，并且使用了诸如自回归综合移动平均值（ARIMA）之类的方法。将在第 4 章，“自回归模型”中进一步讨论。

另一个非常流行的指数平滑是三重指数平滑。三重指数平滑允许您捕获具有水平和趋势的季节性。使用以下等式定义级别，趋势和季节性之间的关系：

在这些等式中，F [t] 捕获了在时间 t 处的观察水平。类似地，T [t] 和 S [t] 在时间 t 捕获趋势和季节性。系数α，β和γ分别表示数据平滑因子，趋势平滑因子和季节性平滑因子，其值在[0,1]之间。这些等式可用于预测下一个时间段，如下所示：

前面的等式可以推广到任何周期 m，如下所示：

术语 F [t] 捕获了季节性成分相对于最近观察到的季节性趋势的偏移。由于将三重指数平滑应用于季节和趋势时间序列，因此可以使用这些值来计算趋势的更好的起始值。让我们使用威斯康星州就业数据演示三阶指数平滑。

可以从这个页面下载威斯康星州就业数据，并可以使用pandas.read_csv命令在 Python 环境中加载 csv 文件。

第一步是加载所需的模块：

# Load modules 
from __future__ import print_function 
import os 
import pandas as pd 
import numpy as np 
from matplotlib import pyplot as plt 

#read the data from into a pandas.DataFrame 
wisc_emp = pd.read_csv('datasets/wisconsin-employment-time-series.csv')

图 3.14：威斯康星州员工数据集

上图显示了威斯康星州员工时间序列数据集。该数据集包括年度趋势和每月季节性。我们知道数据的季节性模式。因此，季节性信息可用于使用以下公式得出趋势的初始值，作为整个季节的平均值：

前面的等式可以在 Python 中实现，如下所示：

# Initialize trend value 
def initialize_T(x, seasonLength): 
    total=0.0 
    for i in range(seasonLength): 
        total+=float(x[i+seasonLength]-x[i])/seasonLength 
    return total

例如，使用以下脚本，由前面的函数生成的初始趋势值是 1.69：

initialize_T(wisc_emp['Employment'], 12)

初始季节性非常关键，可以使用以下函数进行计算：

# Initialize seasonal trend 
def initialize_seasonalilty(x, seasonLength): 
    seasons={} 
    seasonsMean=[] 
    num_season=int(len(x)/seasonLength) 
    # Compute season average 
    for i in range(num_season): 
        seasonsMean.append(sum(x[seasonLength*i:seasonLength*i+seasonLength])/float(seasonLength)) 

    # compute season intial values 
    for i in range(seasonLength): 
        tot=0.0 
        for j in range(num_season): 
            tot+=x[seasonLength*j+i]-seasonsMean[j] 
        seasons[i]=tot/num_season 
    return seasons

将季节的初始值计算为响应 x 的平均值。

一旦获得值，我们就可以设置三重指数预测了：

# Triple Exponential Smoothing Forecast 
def triple_exp_smoothing(x, seasonLength, alpha, beta, gamma, h): 
    yhat=[] 
    S = initialize_seasonalilty(x, seasonLength) 
    for i in range(len(x)+h): 
        if i == 0: 
            F = x[0] 
            T = initialize_T(x, seasonLength) 
            yhat.append(x[0]) 
            continue 
        if i >= len(x): 
            m = i - len(x) + 1 
            yhat.append((F + m*T) + S[i%seasonLength]) 
        else: 
            obsval = x[i] 
            F_last, F= F, alpha*(obsval-S[i%seasonLength]) + (1-alpha)*(F+T) 
            T = beta * (F-F_last) + (1-beta)*T 
            S[i%seasonLength] = gamma*(obsval-F) + (1-gamma)*S[i%seasonLength] 
            yhat.append(F+T+S[i%seasonLength]) 
    return yhat

三重指数平滑由α，β和γ控制。是否存在任何情况都会对结果产生巨大影响。让我们对不同的极端情况进行经验比较，如威斯康星州就业数据集所示：

图 3.15：极限值时的不同配置

使用以下代码生成使用设置为 1 的 alpha，beta 或 gamma 值的三重指数平滑的预测：

# Effect of alpha and beta 
wisc_emp['TEF001'] = triple_exp_smoothing(wisc_emp['Employment'], 12, 0, 0, 1, 0) 
wisc_emp['TEF010'] = triple_exp_smoothing(wisc_emp['Employment'], 12, 0, 1, 0, 0) 
wisc_emp['TEF100'] = triple_exp_smoothing(wisc_emp['Employment'], 12, 1, 0, 0, 0) 

# Plot alpha=0, beta=0, gamma=1 
fig = plt.figure(figsize=(5.5, 5.5)) 
ax = fig.add_subplot(3,1,1) 
wisc_emp['Employment'].plot(color='b', linestyle = '-', ax=ax)  
wisc_emp['TEF001'].plot(color='r', linestyle = '--', ax=ax) 
ax.set_title('TES: alpha=0, beta=0, gamma=1') 

# Plot alpha=0, beta=1, gamma=0 
ax = fig.add_subplot(3,1,2) 
wisc_emp['Employment'].plot(color='b', linestyle = '-', ax=ax)  
wisc_emp['TEF010'].plot(color='r', linestyle = '--', ax=ax) 
ax.set_title('TES: alpha=0, beta=1, gamma=0') 

# Plot alpha=1, beta=0, gamma=0 
ax = fig.add_subplot(3,1,3) 
wisc_emp['Employment'].plot(color='b', linestyle = '-', ax=ax)  
wisc_emp['TEF100'].plot(color='r', linestyle = '--', ax=ax) 
ax.set_title('TES: alpha=1, beta=0, gamma=0') 
fig.subplots_adjust(hspace=.5)

图 3.16：使用（i）alpha：0，beta：0 和 gamma：1 的三重指数平滑进行预测; （ii）alpha：0，beta：1 和 gamma：0；和（iii）alpha：0，beta：0 和 gamma：1

上图显示了使用数据平滑因子 alpha 或季节性因子 gamma 设置为 1 的三重指数平滑。下图显示了使用此脚本针对两个参数设置为一个的相似性绘图：

# Effect of alpha and beta 
wisc_emp['TEF110'] = triple_exp_smoothing(wisc_emp['Employment'], 12, 1, 1, 0, 0) 
wisc_emp['TEF101'] = triple_exp_smoothing(wisc_emp['Employment'], 12, 1, 0, 1, 0) 
wisc_emp['TEF011'] = triple_exp_smoothing(wisc_emp['Employment'], 12, 1, 1, 1, 0)

图 3.17：使用（i）alpha：1，beta：1 和 gamma：0 的三重指数平滑进行预测；（ii）alpha：1，beta：0 和 gamma：1；和（iii）alpha：0，beta：1 和 gamma：1

没有将 beta 设置为零的模型要优于其他模型。让我们使用中间参数运行三次指数平滑，如下图所示：

图 3.18：实际与三次指数平滑比较

可以使用保留样本上的均方根误差来优化预测。让我们比较一下单，双和三重指数平滑的性能，如下图所示：

图 3.19：单指数，双指数和三倍指数平滑之间的比较

根据经验研究，使用平滑或季节性的单一级别可以捕获数据趋势；因此，所有模型均表现出色，因为单指数和双指数能够使用平滑因子进行预测，而三指数平滑能够使用平滑因子或季节性因子来捕获预测。

概括

本章介绍了用于平滑时间序列数据的指数平滑方法。通过包括诸如平滑因子，趋势因子和季节性因子之类的术语，可以轻松地扩展这些方法以进行预测。单阶指数平滑仅使用平滑因子执行平滑，通过包括趋势分量，第二阶平滑因子进一步扩展了平滑因子。还涵盖了三阶平滑处理，该平滑处理将所有平滑处理，趋势和季节性因素合并到模型中。

本章详细介绍了所有这些模型的 Python 实现。平滑方法可用于预测时间序列是否为平稳信号。但是，该假设可能不正确。建议使用高阶指数平滑，但计算起来会很困难。因此，为应对这种方法，提出了其他预测技术，例如 ARIMA，将在下一章中介绍。

四、自回归模型

在上一章中，介绍了基于指数平滑的预测技术，该技术基于以下假设：时间序列由确定性和随机项组成。随机成分为零，并考虑了用于预测的观察数。这假设随机噪声实际上是随机的，并且遵循独立的相同分布。但是，这种假设经常会被违反，平滑化不足以对过程进行建模并建立预测模型。

在这些情况下，自回归模型会非常有用，因为这些模型会利用先前的滞后值通过利用观测值之间固有的序列相关性立即进行调整。本章介绍使用自回归模型的预测概念。自回归模型包括自回归项或移动平均项。根据使用的组件，可以在时间序列预测中使用多种方法，例如移动平均值（MA），自回归移动平均值（ARMA）和自回归综合移动平均值（ARIMA）。本章中的 MA 与第 2 章，“了解时间序列数据”中讨论的移动平均平滑不同。 MA，或更合适的是阶数为q的 MA（q）是基于误差滞后回归的自回归移动平均模型。

本章重点介绍自动回归模型，并将涵盖以下主题：

移动平均线（MA）
自回归（AR）
自回归移动平均线（ARMA）
自回归综合移动平均线（ARIMA）
概括
介绍

时间序列中的自动回归模型的概念是指通过对先前值进行回归而开发的模型。例如，x [t] 是在时间t时的响应，模型的开发如下：

*x[t] = ∅ x[t-1] + Є*

前面的方程式是 AR（1）模型的简单示例。这里，] 是模型系数，Є是误差。此外，类似于回归模型，对于自回归模型以及考虑平稳性或同调误差，误差正态性假设仍然存在。下一个小节介绍了移动平均模型，它是模型与上一个先前值的历史偏差的线性依赖关系。

自回归模型

对时间序列数据进行回归的另一种非常著名的方法是对其滞后项进行回归。这种类型的模型称为自回归模型（AR 模型）。 AR 模型非常适合捕获趋势，因为可以根据先前的时间值预测下一个时间值。因此，AR 模型在下一个预测值是上一个时间段的函数的情况下非常有用，例如由于良好的公司增长而导致的平均股票价格上涨；我们预计随着时间的推移会保留这种影响，价格作为趋势成分随时间的变化将继续增加。

自回归模型定义为 AR（p），其中p是指 AR 组件的顺序。

一阶 AR 模型由 AR（1）表示：

*x[t] = ø∈[t-1] + ∈[t]*

二阶 AR 模型由 AR（2）表示：

*x[t] = ø[1]∈[t-1] + ø[2]∈[t-2] + ∈[t]*

p^第阶 AR 模型由 AR（p）表示：

*x[t] = ø[1]∈[t-1] + ø[2]∈[t-2] + [...] + ø[p]∈[t-p] + ∈[t]*

此处，ø是模型系数，∈ [t] 〜N（0，σ ² ）是时间误差t，并且p是 AR 模型的顺序。让我们利用一种类似于移动平均模型的设置来理解 AR 组件的建模含义。可以使用 statsmodels.tsa 模块中的arma_generate_sample函数生成 AR（1）数据集：

import statsmodels.tsa.api as smtsa 
# Number of samples 
n = 600 
# Generate AR(1) dataset 
ar = np.r_[1, -0.6] 
ma = np.r_[1, 0] 
ar1_data = smtsa.arma_generate_sample(ar=ar, ma=ma, nsample=n)  
plotds(ar1_data)

前面的脚本为 AR（1）场景生成了一个数据集，其中为先前的滞后定义的序列相关性为 0.6。 MA 分量设置为零，以从时间序列信号中消除任何移动平均效果。生成的时间序列信号以及生成的信号的自相关和部分自相关如图 4.1 所示：

图 4.1：AR（1）时间序列信号及其 ACF 和 PACF 图

图 4.7 中的模型关系可以表示为 x [t] =ø [1] x [t-1] +∈ [t] [ 作为数据使用小于 1 的 AR 分量进行模拟；因此，随着ø和t之间的这种关系，自相关将随着时间的推移而降低：

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

因此，ACF 图呈指数下降，而当 PACF 在计算相关性时消除了滞后效应时，仅捕获了有效项。 Φ值会影响信号平稳性。例如，如果在 AR（1）中将Φ从 0.6 增加到 0.95，则模型趋于不稳定，如下图所示：

图 4.2：具有 0.95 的高自相关的 AR（1）时间序列信号

在情形Φ> 1 中，模型变得不稳定。此处显示一个Φ> 1 的非平稳过程的示例：

图 4.3：大于 1 的 AR（1）时间序列信号导致非平稳信号

图 4.9 是在Φ= 1.01 的情况下产生的。类似地，可以生成更高阶的 AR 模型以验证对 PACF 组件的影响。使用以下脚本生成具有 AR（2）和 AR（3）组件的数据集：

# Generate AR(2) dataset 
ar = np.r_[1, 0.6, 0.7] 
ma = np.r_[1, 0] 
ar2_data = smtsa.arma_generate_sample(ar=ar, ma=ma, nsample=n)  
plotds(ar2_data) 

# Generate AR(3) dataset 
ar = np.r_[1, 0.6, 0.7, 0.5] 
ma = np.r_[1, 0] 
ar3_data = smtsa.arma_generate_sample(ar=ar, ma=ma, nsample=n)  
plotds(ar3_data)

图 4.4 和 4.5 显示了时间序列信号及其 AR（2）和 AR（3）时间序列信号的 ACF 和 PACF 图：

图 4.4：AR（2）时间序列及其 ACF 和 PACF

图 4.5：AR（3）时间序列及其 ACF 和 PACF

从不同的数据源可以看出，PACF 正在捕获 AR 分量，q是有效值。可以使用statsmodels.tsa.api模块中的ARMA.fit类评估 AR1 的模型在 Python 中；

# Build AR(1) model 
ar1model = smtsa.ARMA(ar1_data.tolist(), order=(1, 0)) 
ar1=ar1model.fit(maxlag=30, method='mle', trend='nc') 
ar1.summary()

前面的脚本在时间序列数据集上拟合了 AR（1）模型，其结果报告如下：

图 4.6：AR（1）模型输出

使用 0.6 的滞后序列相关性对 AR（1）进行仿真，拟合值的评估值为 0.58，这与实际关系非常接近。同样，将 AR（3）模型拟合到 AR 生成的数据集上，其实际生成的滞后关系分别为 0.6、0.7 和 0.5，以下是拟合值的结果：

图 4.7：AR（3）模型输出

拟合的 AR 模型显示出 0.58、0.67 和 0.44 的关系，这与实际关系非常接近。 AR 和 MA 均可用于校正序列依赖性，但通常，使用 AR 模型可校正正自相关，而使用 MA 模型可校正负相关。

移动平均模型

移动平均模型使用残差之间的依赖关系来预测下一个时间段的值。该模型可帮助您调整任何不可预测的事件，例如导致股市崩溃，导致股价下跌的灾难性事件，这些事件会随着时间的流逝而发生，并被记录为移动平均过程。

MA（1）表示的一阶移动平均值如下：

*x[t] = α - θ[1]Є[t-1] + Є[t]*

MA（2）表示的二阶移动平均值如下：

*x[t] = α - θ[1]Є[t-1] - θ[2]Є[t-2]+ Є[t]*

用 MA（q）表示的q^th 阶移动平均值如下：

*x[t] = α - θ[1]Є[t-1] - θ[2]Є[t-2] - ... - θ[q]Є[t-q]+ Є[t]*

此处，Є [t] 是时间t时的相同独立分布的误差，并且遵循正态分布 N（0，σ ² [Є] ），均值为零，并且σ ² [Є] 方差。 Є [t] 分量表示时间误差t，α和Є标记分别表示平均截距和误差系数。具有q^阶阶的移动平均时间序列模型表示为 MA（q）。前述关系不会更改 MA（q）的期望值，其定义如下：

*E(x[t]) = E(α - θ[1]Є[t-1] - θ[2]Є[t-2] - ... - θ[q]Є[t-q]+ Є[t]) = α*

但是，方差增加，定义如下：

*var(x[t]) = var(α - θ[1]Є[t-1] - θ[2]Є[t-2] - ... - θ[q]Є[t-q]+ Є[t]) = α*

*var(x[t]) = σ<sup class="calibre28">2</sup> (1 + θ[1<sup class="calibre78">2</sup>] + θ[2<sup class="calibre78">2</sup>]+ ... θ[q<sup class="calibre78">2</sup>])*

为了说明移动平均时间序列模型，让我们使用以下代码片段生成信号：

import statsmodels.tsa.api as smtsa 
# Number of samples 
n = 600 
# Generate MA(1) dataset 
ar = np.r_[1, -0] 
ma = np.r_[1, 0.7] 
ma1_data = smtsa.arma_generate_sample(ar=ar, ma=ma, nsample=n)

在前面的代码段中，n表示使用或定义自回归分量的要生成的样本数，ma 解释了该时间的移动平均值分量串联信号。 ar 组件的详细信息将在第 4 章和“自回归模型”的自回归部分中介绍。当前，我们将 ar 对时间序列信号的影响保持为零。上面的代码片段将生成一个时间序列数据集，该数据集具有 MA（1）依赖关系，且误差之间具有 0.7 序列相关性，并且可以表示如下：

*x[t] = 0.7Є[t-1 +] Є[t]*

前面的脚本将生成以下信号：

图 4.8：MA（1）数据集示例

为了评估时间序列信号是由 MA 分量还是 AR 分量组成，使用了自相关（ACF）和部分自相关（PACF）。 ACF 代表：

上式中的分子是时间t和 t-h 的时间序列信号之间的协方差，其中h是时间序列信号中的滞后。除了通过消除间隔之间已经说明的变化来计算相关性外，PACF 的计算方式也与 ACF 类似。这也被定义为条件相关。一阶 ACF 与一阶 PACF 相似。在二阶（滞后）PACF 中，条件概率开始发挥重要作用：

同样，三阶 PACF 可以表示如下：

前面的关系可以进一步扩展到高阶滞后。让我们看一下先前生成的时间序列的 ACF 和 PACF 函数：

图 4.9：MA（1）数据集的 ACF 和 PACF

先前数据集上的 ACF 显示 1 滞后依赖性。由于使用 x [t] =θ∈ [t-1] +∈ [t] 来捕获 MA 关系，因此与滞后项无关， ACF 倾向于捕获 MA 系列的适当顺序q。从图 4.2 中可以看出，ACF 在定义的顺序后不会变为零，而是减小到一个很小的值。置信区间使用关系±2/√n进行测试，其中N和1/√n代表标准偏差的近似值，在独立条件下正确。

使用以下脚本，以q的高阶来查看 ACF 和 PACF 的 MA 组件的影响：

# Generate MA(2) dataset 
ar = np.r_[1, -0] 
ma = np.r_[1, 0.6, 0.7] 
ma2_data = smtsa.arma_generate_sample(ar=ar, ma=ma, nsample=n)  
plotds(ma2_data) 

# Generate MA(3) dataset 
ar = np.r_[1, -0] 
ma = np.r_[1, 0.6, 0.7, 0.5] 
ma3_data = smtsa.arma_generate_sample(ar=ar, ma=ma, nsample=n)  
plotds(ma3_data)

前面的脚本将生成顺序为 MA（2）和 MA（3）的移动平均时间序列信号，而对 AR 分量没有影响。使用上述脚本生成的时间序列信号的 ACF 和 PACF 图如下所示：

图 4.10：MA（2）数据集的示例

图 4.11：MA（3）数据集示例

通常，ACF 很好地定义了错误序列相关性，因此可用于检测 MA（q）；但是，随着阶数的增加以及其他时间序列成分（如季节性，趋势或平稳性）的出现，这使解释变得更加困难。 MA（q）假定过程是平稳的，并且误差是白噪声，以确保无偏估计。例如，在前面的 MA（2）和 MA（3）中，数据中的季节性成分使解释变得更加困难。前面的图是使用以下函数生成的：

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf 
from matplotlib import pyplot as plt 
def plotds(xt, nlag=30, fig_size=(12, 10)): 
    if not isinstance(xt, pd.Series): 
         xt = pd.Series(xt) 
    fig_plt = plt.figure(figsize=fig_size) 
    layout = (2, 2) 

    # Assign axes 
    ax_xt = plt.subplot2grid(layout, (0, 0), colspan=2) 
    ax_acf= plt.subplot2grid(layout, (1, 0)) 
    ax_pacf = plt.subplot2grid(layout, (1, 1)) 

    # Plot graphs 
    xt.plot(ax=ax_xt) 
    ax_xt.set_title('Time Series') 
    plot_acf(xt, lags=50, ax=ax_acf) 
    plot_pacf(xt, lags=50, ax=ax_pacf) 
    plt.tight_layout() 
    return None

可以使用statsmodel.tsa模块中的 ARMA 功能来构建 MA（q）模型。符合 MA（1）模型的示例脚本如下：

# Build MA(1) model 
ma1 = smtsa.ARMA(ma1_data.tolist(), order=(0, 1)).fit( 
    maxlag=30, method='mle', trend='nc') 
ma1.summary()

由于 AR 的阶数保持为零，因此smtsa.ARMA建立了一个 MA（1）。 smtsa.ARMA返回的模型摘要如下所示：

图 4.12：MA（1）模型的输出

如您所见，该模型已捕获残差之间的0.67相关性，这非常接近0.7的模拟值。同样，我们为 MA（3）数据集运行模型，结果如图 4.13 所示：

图 4.13：MA（3）模型在模拟数据集上的输出

使用 ARMA 构建数据集

前两节描述了自回归模型 AR（p），该模型根据自己的滞后项进行回归，移动平均模型 MA（q）建立了误差项的函数。过去的。 AR（p）模型倾向于捕获均值回复效应，而 MA（q）模型倾向于捕获错误的冲击效应，这不是正常事件或不可预测事件。因此，ARMA 模型将 AR 和 MA 组件的功能结合在一起。一个 ARMA（p，q）时间序列预测模型合并了p^th 阶 AR 和 q ^th 分别为 MA 模型。

ARMA（1，1）模型表示如下：

ARMA（1，2）模型表示如下：

ARMA（p，q）模型表示如下：

在此，Φ和θ表示 AR 和 MA 系数。 α和ε [t] 捕获时间 t 的截距和误差。 随着p和q的增加，形式变得非常复杂；因此，滞后算子用于 ARMA 模型的简洁表示。假设L代表滞后算子，并且根据移动的单位，我们将其应用k次。这些运算符也称为后移运算符。

...

使用滞后算子，我们可以如下重写一阶自回归模型：

类似地，运动阶一阶方程可写为：

可以将上述方程式扩展到高阶 AR 和 MA 模型：

可以将前面的两个结合起来形成如下所示的 ARMA：

前面的表示还用于研究冲激响应函数。脉冲响应功能捕获了 x [t] 在时间l（其中 l < t）受到电击时对响应的影响。 在某些外部变化的情况下，脉冲响应也可以被视为对动态系统响应的影响。

让我们通过使用更新的ar和ma组件更新以前使用的脚本来生成 ARMA（1,1）数据集。为了简单起见，我们还将样本数量限制为 600：

# Number of samples 
n = 600 
# Generate AR(1) dataset 
ar = np.r_[1, 0.6] 
ma = np.r_[1, 0.3] 
ar1ma1_data = smtsa.arma_generate_sample(ar=ar, ma=ma, nsample=n)  
plotds(ar1ma1_data )

信号，ACF 和 PACF 图如下所示：

图 4.14：带 ACF 和 PACF 的 ARMA（1,1）信号

通常，销售流程遵循 ARMA（1,1）模型，因为及时销售t是先前发生的及时销售的函数 t-1，在 AR 中发挥作用成分。 ARMA（1,1）的 MA 组件是由公司发起的基于时间的活动引起的，例如，优惠券的分发将导致平均效果在流程中移动，因为销售会暂时增加，并且可以捕捉到销售效果的变化由移动平均线组成。在图 4.14 中，ACF 和 PACF 都显示了正弦曲线，在初始滞后时具有很强的相关性。因此，p和q参数都存在。有多种情况可以选择p和q；可用于确定 ARMA 组件顺序的一些经验法则如下：

自相关呈指数下降，并且 PACF 在滞后 1 处具有显着相关，然后使用p参数
自相关形成正弦波，并且 PACF 在滞后 1 和滞后 2 处具有显着相关，然后对 p 使用二阶值
自相关具有显着的自相关，而 PACF 具有指数衰减，然后存在移动平均，需要设置q参数
自相关显示出显着的串行相关性，而 PACF 显示出正弦波模式，然后设置移动平均值q参数

在 ARMA（1,1）时间序列数据中，由于 ACF 和 PACF 都显示了正弦波模式，p和q这两个参数都影响时间序列信号。对于 ARMA（1,1）时间序列信号，可以使用脉冲响应曲线来计算滞后的影响，如下图所示：

图 4.15：ARMA（1,1）数据集的冲激响应曲线

图 4.15 显示，在五个滞后之后，对响应的影响非常小。为了从数据评估 AR 和 MA 值，使用了statsmodels.tsa.api模块的ARMA.fit函数，如以下脚本所示：

# Build AR(1) model 
ar1ma1 = smtsa.ARMA(ar1ma1_data.tolist(), order=(1, 1)).fit( 
    maxlag=30, method='mle', trend='nc') 
ar1ma1.summary()

图 4.16 显示了上述脚本的输出：

图 4.16：ARMA（1,1）的脉冲响应曲线

模型结果显示 AR 系数值0.58和 MA 值0.29，分别接近 AR 和 MA 组件用来生成时间序列信号的值0.6和0.3。同样，碱性信息标准（AIC）是用于评估模型性能的另一个指标，目的是使 AIC 最小化。要为p和q订单设置数据驱动的评估，可以将 AIC 用作标准。此处显示了 ARMA（1,1）数据集上 AIC 最小化的图示：

# Optimize ARMA parameters
aicVal=[] 
for ari in range(1, 3): 
    for maj in range(1,3): 
        arma_obj = smtsa.ARMA(ar1ma1_data.tolist(), order=(ari, maj)).fit(maxlag=30, method='mle', trend='nc') 
        aicVal.append([ari, maj, arma_obj.aic])

使用p和q作为模型输入，并以 AIC 作为模型输出标准来汇总上述脚本的输出：

表 4.17：ARMA 模型的不同p和q值的 AIC 值

表 4.17显示 ARMA（1,1）是具有最小 AIC 值的最佳模型；但是，ARMA（1,1）和 ARMA（2,2）之间的差异并不高，因为 ARMA（1,1）是自由度（DOF）[HTG1]较小的较简单模型；因此，ARMA（1,1）将比其他复杂模型更受青睐。

让我们说明一个使用实时序列数据的 ARMA 模型。选择用于说明的数据集是 1962 年至 1965 年的 IBM 股票价格数据。第一步是将所需的模块和数据集加载到 Python 环境中：

# Load modules 
from __future__ import print_function 
import os 
import pandas as pd 
import numpy as np 
from matplotlib import pyplot as plt 

# Load Dataset 
ibm_df = pd.read_csv('datasets/ibm-common-stock-closing-prices.csv') 
ibm_df.head() 

#Rename the second column 
ibm_df.rename(columns={'IBM common stock closing prices': 'Close_Price'}, inplace=True) 
ibm_df.head()

前面的脚本利用熊猫来加载数据集。列名使用 pandas DataFrame 支持的rename函数重命名。 IBM 股票收盘价数据集如下所示：

图 4.18：IBM 股票数据集

IBM 图显示了一段时间内数据的显着趋势。该过程的下一步是查看数据集的自相关图。可以使用以下脚本获得 ACF 和 PACF 图：

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf 
ibm_df['residual']=ibm_df['Close_Price']-ibm_df['Close_Price'].mean() 
ibm_df=ibm_df.dropna() 
plot_acf(ibm_df.residual, lags=50) 
plot_pacf(ibm_df.residual, lags=50)

前面的脚本生成的 ACF 和 PACF 图如图 4.16 所示：

图 4.19：来自 IBM 股票数据集的残差的 ACF 和 PACF

ACF 呈线性衰减，显示出很强的序列相关性；然而，部分自相关仅显示了一步依赖性。同样，如图所示，自相关为正，应使用具有一阶相关的 AR 分量进行校正。信号的 QQ 图可帮助您评估正态性假设。 IBM 股票数据集的 QQ 图如图 4.20 所示：

图 4.20：残差的 QQ 图

图 4.21 中的残差图显示了该数据集接近钟形曲线，且拐角案例不遵循正态分布，如图 4.18 所示，为直方图：

图 4.21：IBM 股票收盘价的直方图

为了获得 ARMA 的最佳p和q阶数，使用以下脚本以 AIC 最小化为搜索标准来执行网格搜索：

# Optimize ARMA parameters 
aicVal=[] 
for ari in range(1, 3): 
    for maj in range(0,3): 
        arma_obj = smtsa.ARMA(ibm_df.Close_Price.tolist(), order=(ari, maj)).fit(maxlag=30, method='mle', trend='nc') 
        aicVal.append([ari, maj, arma_obj.aic])

ARMA.fit函数用于使用最大似然准则将 ARMA 预测模型与已定义的p和q模型拟合。模型的结果保存在 aicVal 列表中，如图 4.22 所示：

图 4.22：不同等级 ARMA 模型的 AIC 得分

AIC 建议将 ARMA（1,1）模型作为具有最小 AIC 值的最佳模型。使用以下脚本将 ARMA（1,1）模型重新设置为最佳模型：

# Building optimized model using minimum AIC 
arma_obj_fin = smtsa.ARMA(ibm_df.Close_Price.tolist(), order=(1, 1)).fit(maxlag=30, method='mle', trend='nc') 
ibm_df['ARMA']=arma_obj_fin.predict()

AIC 分数推荐 ARMA（1,0）的模型，AIC 分数为 6702.7。使用 ARMA（1,0）的实际值与拟合值如下：

图 4.23 实际值与预测值

上图显示了 ARMA（1,0）模型的良好拟合，可以预测股票收盘价：

# Plot the curves 
f, axarr = plt.subplots(1, sharex=True) 
f.set_size_inches(5.5, 5.5) 
ibm_df['Close_Price'].iloc[1:].plot(color='b', linestyle = '-', ax=axarr) 
ibm_df['ARMA'].iloc[1:].plot(color='r', linestyle = '--', ax=axarr) 
axarr.set_title('ARMA(1,0)') 
plt.xlabel('Index') 
plt.ylabel('Closing price') Plot the curves 
f, axarr = plt.subplots(1, sharex=True) 
f.set_size_inches(5.5, 5.5) 
ibm_df['Close_Price'].iloc[1:].plot(color='b', linestyle = '-', ax=axarr) 
ibm_df['ARMA'].iloc[1:].plot(color='r', linestyle = '--', ax=axarr) 
axarr.set_title('ARMA(1,1)') 
plt.xlabel('Index') 
plt.ylabel('Closing price')

这些模型的主要限制之一是它们忽略了使信号不稳定的波动性因素。正在考虑的 AR 建模过程是平稳的，即误差项为 IID 且遵循正态分布εt〜N（0，σ ² [ε]）和 |Φ| < 1. |Φ| < 1 条件使时间序列成为有限的时间序列，因为与以前的观测值相比，时间序列中更新的观测值的影响会更高。不满足这些假设的序列属于非平稳序列。此处显示了一个非平稳过程的示例：

图 4.24：AR（1）非平稳过程

从前面的图中可以看出，过程的方差在数据集的末尾处不断增加，并且在 ACF 中观察到了强烈的趋势。下一节讨论的 ARIMA 考虑了对于预测而言非平稳的情况。

有马

ARIMA 也称为 Box-Jenkins 模型，它是 ARMA 模型的概括，它包括集成组件。当数据具有非平稳性时，集成的组件很有用，而 ARIMA 的集成部分有助于减少非平稳性。 ARIMA 对时间序列应用一次或多次差分，以消除非平稳性影响。 ARIMA（ p，d，q ）表示 AR，MA 和差分分量的顺序。 ARMA 模型和 ARIMA 模型之间的主要区别是d组件，该组件更新了构建预测模型的序列。d组件旨在使信号去趋势化以使其稳定，并且 ARMA 模型可以应用于去趋势化的数据集。对于d的不同值，串联响应变化如下：

对于d= 0：x[t] =x[t]

对于d= 1：x[t] =x[t] -x[t-1]

对于d= 2：x[t] =（[x[t] -x[t-1] ）-（x[t-1] -x[t-2] ）=x[t] -2x[t-1] -x[t-2]

从前面的行可以看出，第二个差异不是两个周期之前，而是第一个差异的差异，即d= 1。假设x_hat[t]代表不同的响应，因此 ARIMA 预测可以写成以下形式：

根据p，d和q的顺序，模型的行为会有所不同。例如，ARIMA（1,0，0）是一阶 AR 模型。同样，ARIMA（0,0,1）是一阶 MA 模型。

让我们以 ARIMA（0,1,0）为例说明 ARIMA 建模的不同组成部分。 ARIMA（0,1,0）代表随机游走模型。随机游走模型仅取决于最后一个时间实例，可以表示如下：

*x[t] = x[t-1] + ε[t]*

前面的随机游走方程式也可以用滞后算子表示，如下所示：

*(1-L)x[t] = ε[t]*

在此，εt〜N（0，σ ² ）是误差分量，并且遵循正态分布。在前面的随机游动模型中添加一个常数会导致模型中的漂移，该漂移在本质上也是随机的，如以下等式所示：

*(1-L)x[t] = α + ε[t]*

在此，α是将对时间序列信号产生漂移效果的漂移算子。让我们用 2016 年的道琼斯指数时间序列数据集（DJIA）来说明 ARIMA 建模。DJIA 数据集及其基本 ACF 和 PACF 图如图 4.25 所示：

图 4.25：DJIA 数据集及其 ACF 和 PACF

数据集清楚地显示出具有增加趋势的非平稳信号。 ACF 还显示出指数衰减，而 PACF 在滞后 2 中具有很强的相关性。还可以通过评估不同时间段内的均值和方差来检查非平稳性。均值和方差的差异验证了非平稳性的假设。例如，我们将 DJIA 数据集从 2016 年 1 月至 6 月和 2016 年 7 月至 12 月分成两个学期，并按以下方式评估均值和方差：

mean1, mean2 =djia_df.iloc[:125].Close.mean(), djia_df.iloc[125:].Close.mean() 
var1, var2 = djia_df.iloc[:125].Close.var(), djia_df.iloc[125:].Close.var() 
print('mean1=%f, mean2=%f' % (mean1, mean2)) 
print('variance1=%f, variance2=%f' % (var1, var2)) 

# Output 
mean1=17226.579164, mean2=18616.603593 
variance1=487045.734003, variance2=325183.639530

两个学期的评估均值和方差显示均值和方差值存在显着差异，因此表明数据是不稳定的。评估非平稳性的另一种方法是使用统计测试，例如增强 Dickey-Fuller（ADF）测试。 ADF 是一个单位根测试，用于评估时间序列组件中趋势的强度。 ADF 使用高阶 AR 模型来优化信息标准。假设有一个 AR（3）模型：

可以使用差分滞后项来重写前面的方程式：

让我们将前面的等式简化如下：

可以使用滞后运算符重写前面的公式：

假设 L = 1 是前面多项式的解，则 ADF 求解前面的方程，可以表示如下：

如果 L = 1，则前面的等式可以简化为：

对于单位根，滞后差可以重写为：

使用 Schwartz 贝叶斯信息准则或最小化 Akaike 信息准则（AIC），可以使用上述方程式来确定 AR 滞后分量。在存在强自相关性的情况下，原始序列需要有所区别。 ADF 的 NULL 假设提出 H [0] ：ρ= 0，而对另类假设 H [0] ：ρ< 0。换句话说，零假设是单位根的存在或非平稳性，而替代假设则表明数据的平稳性。

让我们为 DJIA 数据集运行 ADF 测试：

# ADF Test 
from statsmodels.tsa.stattools import adfuller 
adf_result= adfuller(djia_df.Close.tolist()) 
print('ADF Statistic: %f' % adf_result[0]) 
print('p-value: %f' % adf_result[1])

ADF 的输出如下所示：

图 4.26：ADF 测试的输出

理想情况下，ADF 统计信息的负值较大时将代表平稳信号。对于给定的数据集，由于 p 值很高，我们不能拒绝 NULL 假设，使其成为非平稳信号。大多数软件包可确保在执行模型之前满足平稳性。可视化 DJIA 数据集正态性的 qqplot 如图 4.27 所示：

图 4.27：DJIA 数据集的 qqplot

图 4.27 显示了 DJIA 数据集中的显着非正态性。大多数用 Python 编写的核心程序包都会检查数据集中的平稳性。如果不满足平稳性要求，则会引发错误。让我们尝试在当前数据集上强制拟合 ARMA（1,1）模型：

# Force fit ARMA(1,1) model on non-stationary signal 
arma_obj = smtsa.ARMA(djia_df['Close'].tolist(), order=(1, 1)).fit(maxlag=30, method='mle', trend='nc')

先前的力拟合将引发以下错误：

图 4.28：ARMA 模型安装在非平稳信号上时的输出

差异将有助于使信号稳定：

#Let us plot the original time series and first-differences 
first_order_diff = djia_df['Close'].diff(1) 
fig, ax = plt.subplots(2, sharex=True) 
fig.set_size_inches(5.5, 5.5) 
djia_df['Close'].plot(ax=ax[0], color='b') 
ax[0].set_title('Close values of DJIA during Jan 2016-Dec 2016') 
first_order_diff.plot(ax=ax[1], color='r') 
ax[1].set_title('First-order differences of DJIA during Jan 2016-Dec 2016')

图 4.29 显示了上述脚本的输出：

图 4.29：一阶微分后的时间序列信号

图 4.30 显示了d= 1 的积分信号的 ACF 和 PACF 图：

图 4.30：DJIA 数据集的一阶差分 ACF 和 PACF 图

ADF 的残差统计值为-17.13，p-val 接近零，从而表明该模型是平稳的。但是，ACF 和 PACF 都没有显示出移动平均成分呈现随机游走行为的趋势。另外，另一种运行方式是使用 AIC 作为标准进行优化：

# Optimize ARMA parameters 
aicVal=[] 
for d in range(1,3): 
    for ari in range(0, 3): 
        for maj in range(0,3): 
            try: 
                arima_obj = ARIMA(djia_df['Close'].tolist(), order=(ari,d,maj)) 
                arima_obj_fit=arima_obj.fit() 
                aicVal.append([ari, d, maj, arima_obj_fit.aic]) 
            except ValueError: 
                pass

前面的代码以不同的值生成以下 AIC：

图 4.31：不同 ARIMA 模型的 AIC 值

由于模型之间的 AIC 非常接近，因此建议使用主题专家来选择正确的模型。让我们选择 ARIMA（0,2,1）进行模型拟合和评估。 ARIMA（0,2,1）应用二阶微分和一阶移动平均分量来确定观测值之间的关系。设置模型参数，如以下脚本所示：

# Evaluating fit using optimal parameter 
arima_obj = ARIMA(djia_df['Close'].tolist(), order=(0,2,1)) 
arima_obj_fit = arima_obj.fit(disp=0) 
arima_obj_fit.summary() 

# Evaluate prediction 
pred=np.append([0,0],arima_obj_fit.fittedvalues.tolist()) 
djia_df['ARIMA']=pred 
diffval=np.append([0,0], arima_obj_fit.resid+arima_obj_fit.fittedvalues) 
djia_df['diffval']=diffval

使用以下脚本获得与实际值和预测值的比较并将其可视化：

# Plot the curves 
f, axarr = plt.subplots(1, sharex=True) 
f.set_size_inches(5.5, 5.5) 
djia_df['diffval'].iloc[2:].plot(color='b', linestyle = '-', ax=axarr) 
djia_df['ARIMA'].iloc[2:].plot(color='r', linestyle = '--', ax=axarr) 
axarr.set_title('ARIMA(0,2,1)') 
plt.xlabel('Index') 
plt.ylabel('Closing price')

图 4.32：ARIMA（0,2,1）模型的预测输出

图 4.33：使用 ARIMA（0,2,1）模型的 MA 的系数

关于模型是否最优构建的另一个关键检查是白噪声，例如可以使用 qqplot 观察到错误的行为，如图 4.34 所示：

图 4.34：ARIMA（0,2,1）模型的残差的 QQ 正态图

带有 ARIMA（0,2,1）模型的 QQ 正态图显示了显着的正态拟合，如图 4.27 所示。错误正态性也可以使用 Shapiro-wilk 检验进行评估。

ARIMA 模型的进一步扩展包括以大写字母表示的 AR，I 和 MA 的季节性成分。季节性 ARIMA 表示为 ARIMA（p，d，q）（P，D，Q） [m] [，其中 P，D 和 Q 分别代表自回归，积分和移动平均值的季节部分。季节性 ARIMA 模型中的m表示每个季节的周期数。在存在季节性的情况下，可能需要执行额外的季节性差异和季节性调整步骤，以确保信号稳定。例如，如果您查看图 4.30 中显示的 DJIA 差异 ACF 和 PACF 图，则自相关在 42 索引处变得略显重要，这意味着可能存在季节性。季节性差异出现在第一个差异上，可以使用以下代码行看到：

# Seasonality (based on first difference ACF shows significance at 42 lag) 
x=djia_df['Close']-djia_df['Close'].shift(42) 
x.plot()

图 4.35：DJIA 数据集中的季节性

可以使用 statmodel.SARIMAX 模型支持的季节性 ARIMA 来校正先前的季节性。为 DJIA 数据集建立季节性 ARIMA 模型的脚本如下：

# Seasonality (based on first difference ACF shows significance at 42 lag) 
x=djia_df['Close']-djia_df['Close'].shift(42) 
mod = sm.tsa.statespace.SARIMAX(djia_df['Close'], trend='n', order=(0,2,1), seasonal_order=(1,1,1,42)) 
sarimax= mod.fit() 
sarimax.summary()

以下是上述脚本的结果：

图 4.36：SARIMAX 模型的输出

该模型在 AIC 方面显示出显着改进，并且可以针对 SARIMAX 模型中涉及的不同组件进行进一步优化。

置信区间

预测中的常见问题之一是，估计的置信区间是多少？预测模型中的置信度由预测函数中的 alpha 参数定义。 alpha 值 0.05 表示具有 95％置信度的估计，这可以解释为该模型返回的估计具有 5％的概率不落在定义的分布范围内。置信区间的评估如下：

在此，Z [α] 是基于 alpha 定义的临界值。对于 alpha 值 0.05，临界值为 1.96。可以使用arima_obj_fit对象的预测函数，获得使用 ARIMA（0,2,1）模型建模的 DJIA 数据集的 alpha 值为 0.05 的置信区间：

# Forecasting and CI 
f, err, ci=arima_obj_fit.forecast(40) 
plt.plot(f) 
plt.plot(ci) 
plt.xlabel('Forecasting Index') 
plt.ylabel('Forecasted value')

使用上述脚本获得的预测估计值和置信区间如图 4.37 所示：

图 4.37：95％置信度的预测和置信区间

概括

在本章中，我们介绍了自回归模型，例如 MA 模型，以利用误差关系捕获序列相关性。在类似的直线上，涵盖了 AR 模型，该模型使用滞后作为相关观测值来建立预测。 AR 模型非常适合捕获趋势信息。还说明了基于 ARMA 的方法，该方法集成了 AR 和 MA 模型以捕获任何基于时间的趋势和灾难性事件，从而导致许多错误，而这些错误将需要时间进行纠正，例如经济崩溃。所有这些模型都假设是平稳的。在不存在平稳性的情况下，提出了基于差异的模型（例如 ARIMA），该模型在时间序列数据集中执行差异处理以删除任何与趋势相关的成分。通过使用 Python 的 tsa 模块的示例说明了预测方法。

本章重点介绍使用统计方法进行预测。下一章将统计方法扩展到机器学习方法，以专门针对深度学习模型进行预测。

五、用于进行时间序列预测的深度学习

到目前为止，在本书中，我们已经描述了用于时间序列分析的传统统计方法。在前面的章节中，我们讨论了从过去的观察结果中预测未来某个时间序列的几种方法。一种进行预测的方法是自回归（AR）模型，该模型将时间t的序列表示为先前的线性回归 p 观察结果：

Є [t] 是 AR 模型的残差项。

可以将线性模型的基本思想概括为：时间序列预测的目的是开发一个函数f，该函数根据 x [t] 的预测值进行预测先前的p时间点：

*x[t] = f(x[t-1],x[t-2], ... ,x[t-p])*

在本章中，我们将探索三种基于神经网络的方法来开发函数f。每种方法都包括定义一个神经网络体系结构（根据隐藏层的数量，每个隐藏层中神经元的数量等等），然后使用反向传播算法或其适合于当前网络体系结构的变体来训练网络。用过的。

最近几年见证了对神经网络的兴趣兴起。由于可以从数字媒体获得大量训练数据，并且可以更便宜地访问 GPU 驱动的并行计算，因此这之所以成为可能。这些因素使训练神经网络具有数十万个参数，在某些情况下甚至具有数百万个参数。不同体系结构的神经网络已成功应用于解决计算机视觉，语音识别和自然语言翻译方面的问题。在这些领域中设计和训练神经网络的研究和实践被普遍称为深度学习，顾名思义，它指示了这些模型中使用的许多隐藏层。

深度学习提供了有趣的神经体系结构，旨在解决图像和语言数据的特殊结构特征。例如，卷积神经网络（CNN）旨在利用图像的二维或三维结构，而大多数语言模型都使用递归神经网络（RNN）支持口语和书面语言固有的顺序和记忆。这些新的发展也已应用于传统上统计机器学习占主导地位的领域。这样的领域之一是时间序列预测。

在本章中，我们将探讨三种不同类型的神经网络，用于时间序列预测。我们从 Multi-layer perceptrons 开始（mlp）。这将是经常性的神经网络，其适用于数据点的连续排列。最后，我们将涵盖卷积神经网络，这些网络主要用于图像，但我们将讨论 CNN 的特殊形式如何用于时间序列预测。这些主题将通过对网络架构的解释以及如何应用于时间序列预测来涵盖这些主题。代码演示向您展示如何使用尖端深度学习库来开发时间序列预测模型。本章中的示例是使用 Keras API 来实现的，用于深入学习。 Keras 是一个高级 API，允许定义不同的神经网络架构并使用各种基于梯度的优化器训练它们。在后端，Keras 使用在 C，C ++ 和 Fortran 中实现的低级计算框架。有几种这样的低级框架是可用的开源。 Keras 支持以下三个：Tensorflow 由 Google 开发，是 Keras，CNTK 的默认后端，来自 Microsoft，以及 Theano 的开源框架，最初在加拿大蒙特利尔大学开发。本书中的示例使用 TensorFlow 作为后端。因此，要运行示例，您需要安装 keras 和 tensorflow。

因此，不用费劲，让我们深入研究深度学习中的时间序列预测的迷人主题。

多层感知器

多层感知器（MLP）是神经网络的最基本形式。 MLP 由三个组件组成：输入层，一堆隐藏层和输出层。输入层表示回归变量或输入特征的向量，例如，从先前的p时间点[ x [t-1] ，x [t-2] ，...，x [tp] ]。输入要素被馈送到具有n神经元的隐藏层，每个神经元对输入要素应用线性变换和非线性激活。神经元的输出为 g [i] = h（ w [i] x + b [i] ），其中 w [i] b [i] 是线性变换和 h [ 是非线性激活函数。非线性激活函数使神经网络能够对回归变量与目标变量之间的基础关系的复杂非线性进行建模。通常，h是 S 型函数1 / (1 - exp(-z))，它将任意实数压缩为区间 [0,1]。由于具有此特性，因此会使用 S 型函数来生成二进制类别的概率，因此通常在分类模型中使用。非线性激活函数的另一种选择是tanh函数(1 - exp(-z)) / (1 + exp(-z)))，它将任何实数绑定到区间[-1, 1]。在某些情况下，h是恒等式或线性函数。

在单个隐藏层神经网络的情况下（如下图左侧所示），每个神经元的输出将传递到输出层，该输出层将应用线性变换和激活函数来生成目标变量的预测在时间序列预测的情况下，是在 t ^th 时间点的序列的预测值。在 MLP 中，多个隐藏层相互堆叠。来自一个隐藏层的神经元输出将作为输入馈送到下一个隐藏层。该隐藏层中的神经元会转换输入并传递到下一个隐藏层。最后，最后一个隐藏层进入输出层：

图 5.1：多层感知器

MLP 的隐藏层也称为密集或有时完全连接的层。 “密集”一词具有以下含义：密集层的所有神经元都与上一层和下一层中的所有神经元相连。如果前一层是输入层，则所有输入要素都将馈入隐藏层的每个神经元。由于输入层和第一密集层之间以及密集层本身之间的多对多连接，MLP 具有大量可训练的权重。例如，如果输入特征的数量为p，并且存在具有神经元数量 n [1]，n [2 的三个致密层] 和 n [3]，则可训练砝码的数量为p × n[1] + n[1] × n[2] + n[2] × n[3] + n[3]。此计算中的最后一个元素是连接第三隐藏层和输出层的权重数。深度 MLP 具有几层密集层，每层中有数百甚至数千个神经元。因此，深层 MLP 中可训练权重的数量非常大。

培训 MLP

通过运行基于梯度的优化算法（例如随机梯度下降）来找到神经网络的权重w，该算法可迭代地使网络在进行预测时产生的损失或误差（L）最小化在训练数据上。 均方误差（MSE）和平均绝对误差（MAE）（有时表示绝对百分比误差）通常用于回归二进制和分类日志损失是分类问题的常见损失函数。对于时间序列预测，MSE 和 MAE 将易于训练神经模型。

梯度下降算法通过沿权重的梯度路径移动权重i来工作。梯度是损失函数L相对于重量的偏导数。最简单的更改权重w的更新规则需要权重的值，L相对于权重的偏导数以及控制点下降速度的学习率α 沿渐变：

此基本更新规则具有多种变体，会影响算法的收敛性。但是，所有基于梯度的算法的关键输入是必须为网络的所有权重计算的偏导数。在深度神经网络（其中一些具有数百万的权重）中，导数计算可能是庞然大物的计算任务。这正是著名的反向传播算法有效解决这一问题的地方。

要了解反向传播，首先应该了解计算图以及它们如何用于在神经网络中进行计算。

让我们考虑一个简单的单隐藏层神经网络，该网络具有两个隐藏单元，每个隐藏单元均具有 S 型激活。输出单元是其输入的线性变换。网络中有两个输入变量[ x [1] ，x [2] ]。权重沿网络的边缘显示：

网络执行一系列加法，乘法和几个 S 形函数，以将输入转换为预测y_hat。通过神经网络将输入转换为预测的过程称为前向通过。下图显示了如何通过计算图对输入对 [-1,2] 实现正向传递。每次计算都会产生一个中间输出 p [ i ]。中间结果 p [ 7 ] 和 p [ 8 ] 是隐藏神经元 g [ 1 ] [ 和 g [ 2 ]。在训练过程中，损耗L也在前向通过中计算。

图 5.3：具有两个隐藏神经元的单层感知器的计算图

此时，将采用反向传播算法来计算通过边连接的两个节点之间的偏导数。图中用于计算偏导数的向后遍历也称为向后遍历。将偏微分算子应用于每个节点，并将偏导数分配给沿计算图连接下游节点的各个边。遵循链规则，通过将连接权重节点和损失节点的所有边上的偏导数相乘来计算偏导数∂L/∂w。如果权重节点和损失节点之间存在多条路径，则将沿每条路径的偏导数相加，以获得损失相对于权重的总偏导数。这种基于图形的实现正向和向后传递的技术是强大的深度学习库中使用的基础计算技巧。下图显示了反向传递：

图 5.4：计算图中的偏导数计算

损失相对于权重的偏导数可通过应用链式法则获得：

在训练期间，权重使用通常从具有[ -1,1 ]的上限和下限的均匀分布或均值为零和单位方差的正态分布采样的随机数初始化。这种随机初始化方案具有一些增强优化收敛性的变体。在这种情况下，我们假设权重是根据均匀随机分布进行初始化的，因此 w [1] = -0.33，w [2] = 0.57，w [3] = 0.02，w [4] = -0.01，w [5] = 0.07 和 w [6] = 0.82。使用这些值，让我们遍历正向和向后遍历计算图。我们用蓝色的前向通过期间计算的值和红色的反向通过期间计算的梯度来更新前一图。在此示例中，我们将目标变量的实际值设置为 y = 1：

图 5.5：向前（蓝色）和向后（以红色）通过计算图

一旦计算了沿边缘的梯度，相对于权重的偏导数就只是链式规则的应用，我们之前已经讨论过。偏导数的最终值如下：

下一步是使用梯度下降算法更新权重。因此，在学习率为α = 0.01 的情况下，w [5] 的新值= 0.07-0.01 x -0.384 = 0.0738。其余权重也可以使用类似的更新规则进行更新。

迭代权重更新的过程重复多次。权重更新完成的次数称为训练数据的经过次数或通过次数。通常，与前一时期相比，损失函数变化的公差标准控制着时期的数量。

反向传播算法与基于梯度的优化器一起用于确定神经网络的权重。值得庆幸的是，存在强大的深度学习库，例如 Tensorflow，Theano 和 CNTK，它们实现了计算图来训练任何体系结构和复杂性的神经网络。这些库具有内置支持，可以将计算作为对多维数组的数学运算来运行，并且还可以利用 GPU 进行更快的计算。

用于时间序列预测的 MLP

在本节中，我们将使用 MLP 来开发时间序列预测模型。这些示例中使用的数据集是关于空气污染的，该污染是通过直径小于或等于 2.5 微米的颗粒物（PM）的浓度测量的。还有其他变量，例如气压，空气温度，露点等。已经开发了几个时间序列模型，一个是气压，另一个是 pm 2.5。该数据集已从 UCI 机器学习存储库下载。问题描述和数据集的链接为这个页面。

气压时间序列模型的代码在 Jupyter 笔记本中code/Chapter_5_Air Pressure_Time_Series_Forecasting_by_MLP.ipynb，而pm2.5上的代码在code/ Chapter_5_Air Pressure_Time_Series_Forecasting_by_MLP.ipynb中。 code文件夹位于该书的 GitHub 存储库中。现在让我们描述如何开发气压的时间序列模型。

我们首先导入运行代码所需的软件包：

from __future__ import print_function 
import os 
import sys 
import pandas as pd 
import numpy as np 
%matplotlib inline 
from matplotlib import pyplot as plt 
import seaborn as sns 
import datetime

设置当前工作目录后，将数据集从.csv文件读取到pandas.DataFrame：

#set current working directory 
os.chdir('D:/Practical Time Series') 
#Read the dataset into a pandas.DataFrame 
df = pd.read_csv('datasets/PRSA_data_2010.1.1-2014.12.31.csv')

为确保行按观察日期和时间的正确顺序，从与日期和时间相关的列中创建一个新列datetime。新列由 Python 的datetime.datetime对象组成。 DataFrame在此列上按升序排序：

df['datetime'] = df[['year', 'month', 'day', 'hour']] 
                 .apply(lambda row: datetime.datetime(year=row['year'], month=row['month'], day=row['day'], 
                                                                                          hour=row['hour']), axis=1) 
df.sort_values('datetime', ascending=True, inplace=True)

PRES列包含有关气压的数据。让我们绘制一个箱形图以可视化PRES的集中趋势和离散度：

plt.figure(figsize=(5.5, 5.5)) 
g = sns.boxplot(df['PRES']) 
g.set_title('Box plot of Air Pressure')

尽管箱形图未显示实际的时间序列，但它对于快速识别异常值很有用。按照惯例，位于第 25 ^个四分位数之外的观测值是四分位数间距离的 1.5 倍，而在第^个四分位数的范围内是 75 四分位数的范围外+ 1.5 倍。四分位数之间的范围是第 75 ^个四分位数–第 25 ^个四分位数。异常值的存在可用于确定训练神经网络的损失函数，正如我们将在本章的示例中看到的那样。

气压的箱线图显示没有观测值可以指定为离群值。

图 5.6：气压箱线图

接下来，我们绘制整个观察期间的气压时间序列：

g = sns.tsplot(df['PRES']) 
g.set_title('Time series of Air Pressure') 
g.set_xlabel('Index') 
g.set_ylabel('Air Pressure readings in hPa')

图 5.7：气压的时间序列

如果变量在[ -1,1 ]范围内，则梯度下降算法的效果更好（例如，收敛速度更快）。许多源将边界放宽到[ -3,3 ]。将PRES变量进行最小最大缩放，以将转换后的变量限制在[0,1]内。变量x的最小值最大值缩放按以下方式进行：

前面的公式将导致x_scaled∈[0,1]：

from sklearn.preprocessing import MinMaxScaler
scaler = MinMaxScaler(feature_range=(0, 1))
df['scaled_PRES'] = scaler.fit_transform(np.array(df['PRES'])
 .reshape(-1, 1))

变量scaled_PRES用于开发时间序列预测模型。为了对原始变量PRES进行预测，我们对模型的预测应用了适当的逆变换。

在培训模型之前，数据集分为两部分，列车集和验证集。神经网络在火车集上培训。这意味着在列车集上完成损耗函数，反向衰减和权重更新的计算。验证集用于评估模型并确定培训的时期数量。增加时期的数量将进一步降低火车集的损失函数，但由于火车集上的过度装备而可能不一定对验证集具有相同的效果。因此，通过对验证集计算的损耗函数保持挖掘来控制时代的数量。我们使用带有 TensoRFlow 后端的 Keras 来定义和培训模型。模型训练和验证中涉及的所有步骤都是通过调用 Keras API 的适当函数来完成的。

首先，将数据分为训练集和验证集。数据集的时间段为

2010 年 1 月 1 日至 2014 年 12 月 31 日。前四年（2010 年至 2013 年）用作火车和

保留 2014 以供验证：

split_date = datetime.datetime(year=2014, month=1, day=1, hour=0) 
df_train = df.loc[df['datetime']<split_date] 
df_val = df.loc[df['datetime']>=split_date] 
print('Shape of train:', df_train.shape) 
print('Shape of test:', df_val.shape)

前面的代码行生成以下输出：

Shape of train: (35064, 15) 
Shape of test: (8760, 15)

还绘制了scaled_PRES的训练和验证时间序列：

plt.figure(figsize=(5.5, 5.5))
g = sns.tsplot(df_train['scaled_PRES'], color='b')
g.set_title('Time series of scaled Air Pressure in train set')
g.set_xlabel('Index')
g.set_ylabel(‘Scaled Air Pressure readings')

图 5.8：训练数据中最小最大比例气压的时间序列

plt.figure(figsize=(5.5, 5.5)) 
g = sns.tsplot(df_val['scaled_PRES'], color='r') 
g.set_title('Time series of standardized Air Pressure in validation set') 
g.set_xlabel('Index') 
g.set_ylabel('Standardized Air Pressure readings')

图 5.9：验证数据中最小最大标定气压的时间序列

现在我们需要生成回归变量（X）和目标变量（y），以进行训练和验证。从DataFrames中scaled_PRES列的原始 1-D 数组创建一个 2D 回归数组和目标的 1-D 数组。对于时间序列预测模型，过去 7 天的观测值用于预测第二天。这等效于 AR（7）模型。我们定义了一个函数，该函数以原始时间序列和时间步长为输入，以返回X和y的数组：

def makeXy(ts, nb_timesteps): 
    """ 
    Input:  
           ts: original time series 
           nb_timesteps: number of time steps in the regressors 
    Output:  
           X: 2-D array of regressors 
           y: 1-D array of target  
    """ 
    X = [] 
    y = [] 
    for i in range(nb_timesteps, ts.shape[0]): 
        X.append(list(ts.loc[i-nb_timesteps:i-1])) 
        y.append(ts.loc[i]) 
    X, y = np.array(X), np.array(y) 
    return X, y

调用makeXy函数来生成训练和验证集：

X_train, y_train = makeXy(df_train['scaled_PRES'], 7) 
print('Shape of train arrays:', X_train.shape, y_train.shape) 
X_val, y_val = makeXy(df_val['scaled_PRES'], 7) 
print('Shape of validation arrays:', X_val.shape, y_val.shape)

先前的print函数的输出如下：

Shape of train arrays: (35057, 7) (35057,) 
Shape of validation arrays: (8753, 7) (8753,)

现在，我们使用 Keras 功能 API 定义 MLP。在这种方法中，各层被声明并级联为彼此的输入和输出：

from keras.layers import Dense, Input, Dropout 
from keras.optimizers import SGD 
from keras.models import Model 
from keras.models import load_model 
from keras.callbacks import ModelCheckpoint

输入层以形状(None, 7)和类型float32声明。 None指示在运行时确定的实例数：

input_layer = Input(shape=(7,), dtype='float32')

密集层通过线性激活来声明：

dense1 = Dense(32, activation='linear')(input_layer) 
dense2 = Dense(16, activation='linear')(dense1) 
dense3 = Dense(16, activation='linear')(dense2)

多个隐藏层和每个隐藏层中的大量神经元使神经网络能够对回归变量与目标之间基本关系的复杂非线性进行建模。但是，深度神经网络也可能过度拟合训练数据，并在验证或测试集上给出较差的结果。辍学已被用来规范化深度神经网络。在此示例中，在输出层之前添加了一个辍学层。删除操作会在进入下一层之前，将输入神经元的 p 分数随机设置为零。随机丢弃的输入本质上充当了模型集合的自举聚合或袋装类型。随机森林通过在输入要素的随机子集上构建树来使用装袋。我们使用p = 0.2删除 20％的随机选择输入要素：

dropout_layer = Dropout(0.2)(dense3)

最后，输出层可以预测第二天的气压：

output_layer = Dense(1, activation='linear')(dropout_layer)

输入，密集和输出层现在将打包在Model中，这是用于训练和做出预测的包装器类。 均方误差（MSE）被用作loss功能。

网络权重通过 Adam 算法优化。亚当代表自适应矩估计，已经成为训练深度神经网络的流行选择。与随机梯度下降不同，Adam 对每个权重使用不同的学习率，并随着训练的进行分别更新。权重的学习率基于权重梯度和平方梯度的指数加权移动平均值进行更新：

ts_model = Model(inputs=input_layer, outputs=output_layer) 
ts_model.compile(loss='mean_squared_error', optimizer='adam') 
ts_model.summary()

summary功能显示分层详细信息，例如输入和输出的形状以及可训练砝码的数量：

Layer (type)                 Output Shape              Param #   =================================================================
input_6 (InputLayer)         (None, 7)                 0         _________________________________________________________________
dense_21 (Dense)             (None, 32)                256       _________________________________________________________________
dense_22 (Dense)             (None, 16)                528       
_________________________________________________________________
dense_23 (Dense)             (None, 16)                272       _________________________________________________________________
dropout_6 (Dropout)          (None, 16)                0         _________________________________________________________________
dense_24 (Dense)             (None, 1)                 17        
=================================================================
Total params: 1,073
Trainable params: 1,073
Non-trainable params: 0

通过在模型对象上调用fit函数并传递**X_train**和**y_train**来训练模型。训练针对预定数量的时期完成。另外，**batch_size**定义了用于反向传播实例的火车样本集数量。在每个时期完成后，还将传递验证数据集以评估模型。 **ModelCheckpoint**对象跟踪验证集上的损失函数，并为损失函数最小的时期保存模型：

save_weights_at = os.path.join('keras_models', 'PRSA_data_Air_Pressure_MLP_weights.{epoch:02d}-{val_loss:.4f}.hdf5') 
save_best = ModelCheckpoint(save_weights_at, monitor='val_loss', verbose=0, 
                            save_best_only=True, save_weights_only=False, mode='min', 
                            period=1) 
ts_model.fit(x=X_train, y=y_train, batch_size=16, epochs=20, 
             verbose=1, callbacks=[save_best], validation_data=(X_val, y_val), 
             shuffle=True)

Jupyter 笔记本code/ Chapter_5_Air Pressure_Time_Series_Forecasting_by_MLP.ipynb中详细介绍了训练集和验证集上的按时代划分的 MSE。还显示了完成每个时期所需的时间。

预测是根据保存得最好的模型进行的。对模型中基于标准气压的预测进行逆变换，以获得对原始气压的预测。拟合优度或 R 平方也被计算：

best_model = load_model(os.path.join('keras_models', 'PRSA_data_Air_Pressure_MLP_weights.06-0.0039.hdf5')) 
preds = best_model.predict(X_val) 
pred_PRES = mu + sigma*preds 
pred_PRES = np.squeeze(pred_PRES) 

from sklearn.metrics import r2_score 
r2 = r2_score(df_val['PRES'].loc[7:], pred_PRES) 
print('R-squared for the validation set:', round(r2,4))

验证集上模型的 R 平方为 0.9957。最后，绘制了前五十个实际和预测的气压值：

plt.figure(figsize=(5.5, 5.5)) 
plt.plot(range(50), df_val['PRES'].loc[7:56], linestyle='-', marker='*', color='r') 
plt.plot(range(50), pred_PRES[:50], linestyle='-', marker='.', color='b') 
plt.legend(['Actual','Predicted'], loc=2) 
plt.title('Actual vs Predicted Air Pressure') 
plt.ylabel('Air Pressure') 
plt.xlabel('Index')

图 5.10：实际和 MLP 预测的气压时间序列

还为pm2.5变量训练了使用 MLP 的时间序列预测模型，详细实现在 Jupyter 笔记本code\Chapter_5_PM2.5_Time_Series_Forecasting_by_MLP.ipynb中。

绘制了pm2.5的箱线图，以检查异常值的存在：

plt.figure(figsize=(5.5, 5.5)) 
g = sns.boxplot(df['pm2.5']) 
g.set_title('Box plot of pm2.5')

图 5.11：PM2.5 时间序列的箱线图

如上图所示，pm2.5具有异常值，因此，选择 MSE 作为训练 MLP 的损失函数并不是最佳方法。 MSE 是实际值与预测值之间的偏差的平方，它会给损失函数带来巨大的波动。这使梯度下降算法不稳定，并对其收敛产生不利影响。 MAE 绝对值是一阶差，因此不太容易受异常值的影响。因此，在这种情况下，MAE 用于训练神经网络。

在继续定义和训练 MLP 之前，让我们仔细查看一年多和六个月以上的时间序列，以查看是否明显存在任何模式，例如趋势，季节性等。

plt.figure(figsize=(5.5, 5.5)) 
g = sns.tsplot(df['pm2.5'].loc[df['datetime']<=datetime.datetime(year=2010,   month=6,day=30)], color='g') 
g.set_title('pm2.5 during 2010') 
g.set_xlabel('Index')  
g.set_ylabel('pm2.5 readings')

plt.figure(figsize=(5.5, 5.5))
g = sns.tsplot(df['pm2.5'].loc[df['datetime']<=datetime.datetime(year=2010, month=1,day=31)], color='g')
g.set_title('pm2.5 during Jan 2010')
g.set_xlabel('Index')
g.set_ylabel('pm2.5 readings')

图 5.12：2010 年 1 月至 2010 年 6 月 PM2.5 的时间序列

图 5.13：2010 年 1 月至 2010 年 12 月期间 PM2.5 的时间序列

当我们放大pm2.5的六个月数据时，就会看到一些模式。 pm2.5的时间序列显示出周期性的波峰和波谷，尽管两个高点和两个低点之间的时间间隔有所不同。波峰的高度显示出很大的波动。此外，波峰和波谷都分布在多个时间步长上。这由被几个较小的峰包围的高峰表示。同样，槽中的波动也很小。如下图所示，整个系列（从 2010 年到 2014 年）没有显示任何长期趋势，尽管整个系列中可能都存在短期趋势：

plt.figure(figsize=(5.5, 5.5))
g = sns.tsplot(df['pm2.5'])
g.set_title('Time series of pm2.5')
g.set_xlabel('Index')
g.set_ylabel('pm2.5 readings')

图 5.14：2010 年 1 月至 2014 年 12 月期间 PM2.5 的时间序列

具有多个隐藏层和每个隐藏层中多个神经元的神经网络将适合于对数据中的这种复杂非线性模式进行建模。让我们尝试使用 MLP。

**pm2.5**的 MLP 分别具有三层致密层，分别在第一，第二和第三层中具有三十二，十六和十六个隐形神经元。每层都有tanh激活。输出层具有线性激活。使用亚当优化器与 MAE 作为损失功能进行培训：

input_layer = Input(shape=(7,), dtype='float32') 

dense1 = Dense(32, activation='tanh')(input_layer) 
dense2 = Dense(16, activation='tanh')(dense1) 
dense3 = Dense(16, activation='tanh')(dense2) 

dropout_layer = Dropout(0.2)(dense3) 

output_layer = Dense(1, activation='linear')(dropout_layer) 

ts_model = Model(inputs=input_layer, outputs=output_layer) 
ts_model.compile(loss='mean_absolute_error', optimizer='adam')

最佳模型的验证集给出的 MAE 为 11.8993。下图绘制了前五十个实际值和预测值：

图 5.15：PM2.5 的实际时间和 MLP 预测的时间序列

递归神经网络

到目前为止，我们已经使用 MLP 来开发时间序列预测模型。为了预测时间[ x [t-1] ，...，x [tp] ]时的序列p，我们提供了过去p的输入向量时间步长p到 MLP。过去的p时间步长作为不相关的独立变量输入到 MLP。这种模型的一个问题是，它没有隐式考虑观察值彼此相关的时间序列数据的顺序性质。时间序列中的相关性也可以解释为时间序列自身的记忆。在本节中，我们将讨论递归神经网络（RNN），它们在结构上与 MLP 不同，更适合于顺序数据。

RNN 成为开发语言模型的好选择，这种语言模型可以根据出现在单词之前的单词来模拟单词出现的概率。到目前为止，RNN 已用于开发进行文本分类的模型，例如情感预测，语言翻译，并与卷积神经网络结合使用，通过文本生成来描述图像。

下图显示了具有x_hat[t]时间步长的 RNN，每个时间步长都由在相应步长处出现的输入提供。最后一个时间步的输出是输入序列的预测。例如，此 RNN 可用于开发时间序列预测模型，其中输入序列[ x [t-1] ，...，x [ tp ] ] 被馈送到 RNN，最后一个时间步的输出是预测：

图 5.16：具有p时间步长的递归神经网络

RNN 的每个计算单元的内部状态为s[0]，并承载该系列的内存。每个时间步长中隐藏神经元的数量是内部状态的维数，其计算方式为f。每个时间步长的这种内部状态都承载着该系列的记忆。对于第一个时间步，输入y = h(V s[t-1])初始化为全零。函数h(·)是非线性变换，例如tanh。最后一个时间步长的输出是最后一个内部状态U的函数，其中W是非线性变换。还可以使 RNN 返回每个时间步的输出，这对于要求翻译后的文本为目标语言中的一系列单词的语言翻译模型很有用。要注意的重要一点是，权重V，x[t]和x[t] RNN 的时间步中共享，因此网络的可训练参数数量保持较低。

在语言模型的情况下，输入m可以是单词的一键编码。对于一个变量的时间序列建模，x[t]是一个数字。但是，RNN 可以应用于多元时间序列，这对于解决各个序列之间的互相关性特别有用。如果 RNN 对m时间序列建模，则 x [ i ] 是m维向量。

双向递归神经网络

到目前为止讨论的 RNN 是单向的，并且沿原始时间序列的方向遍历。但是，在许多情况下，捕获反向序列中的顺序信息和内存可以改善预测。这种同时使用前向和后向遍历的 RNN 称为双向 RNN，它可以提高网络在长距离内捕获内存的能力。下图显示了双向 RNN：

图 5.17：具有p时间步长的双向递归神经网络

深度递归神经网络

深度学习的强大功能是将多个计算层彼此堆叠在一起。对于 MLP，多个隐藏层彼此相对放置。我们可以通过将多个 RNN 相互堆叠来制作深层 RNN。在深度 RNN 中，递归层的输入序列是前一个递归层的输出序列。最终预测是从最终 RNN 层的最后时间步获得的。下图说明了深层的 RNN：

图 5.18：具有p时间步长的深度递归神经网络

我们也可以通过将多个双向 RNN 相互堆叠来创建深度双向 RNN。这种深层的 RNN 用于复杂的任务，例如语言翻译和文本生成，以描述可以是图像的上下文。显然，深度 RNN 还需要数百万的网络权重才能得到训练。

训练递归神经网络

众所周知，RNN 很难训练。到目前为止，我们一直在谈论的 RNN 属于 Vanilla RNN，其梯度逐渐消失和爆炸，导致训练期间的结果不稳定。结果，RNN 难以学习远程依赖关系。对于时间序列预测，过去追溯太多的时间步将是有问题的。为了解决此问题，长短期记忆（LSTM）和门控循环单元（GRU）是 RNN 的特殊类型。，已介绍。在本章中，我们将使用 LSTM 和 GRU 开发时间序列预测模型。在此之前，让我们回顾一下如何使用时间反向传播（BPTT）训练 RNN，这是反向传播算法的一种变体。我们将发现在 BPTT 期间消失梯度和爆炸梯度是如何产生的。

让我们考虑一下我们用于时间序列预测的 RNN 的计算图。梯度计算如下图所示：

图 5.19：具有 p 个时间步长的深度递归神经网络的时间反向传播

对于权重U，有一条路径可以计算偏导数，如下所示：

但是，由于 RNN 的顺序结构，存在多个连接权重和损失以及与损失的路径。因此，偏导数是沿着各个路径的偏导数的总和，该总导数从损耗节点开始，到计算图中的每个时间步节点结束：

通过对连接损失节点和每个时间步节点的路径求和来计算权重梯度的技术是沿时间反向传播，这是原始反向传播算法的特例。远程 RNN 中梯度消失的问题是由于 BPTT 梯度计算中的乘法项引起的。现在，让我们从上述方程式中的一个检查乘法项。

沿着连接损耗节点和第^个第个时间步长的计算路径的梯度为：

梯度乘法的链很长，很难建模远距离依存关系，这就是梯度消失的问题所在。

内部状态[0, 1)的激活功能为tanh或 S 型。 tanh的一阶导数是exp(-x) / (1 + exp(-x))²，它绑定在(0, 1/4]中。对于 S 型函数，一阶导数为∂s[i]/∂s[i-1]，它绑定在s[i]中。因此，梯度W为正分数。对于长时程，将这些分数梯度相乘会使最终乘积减小为零，并且长时步中没有梯度流。由于梯度的值可以忽略不计，权重不会更新，因此据说神经元已饱和。

值得注意的是，U，∂s[i]/∂s[i-1]和∂s[i]/∂W是矩阵，因此偏导数t和s[t]是在矩阵上计算的。最终输出通过矩阵乘法和加法计算。矩阵的一阶导数称为 Jacobian。如果雅可比矩阵的任何一个元素都是分数，那么对于远程 RNN，我们将看到消失的梯度。另一方面，如果雅可比行列式的元素大于 1，则训练过程会遭受爆炸梯度的影响。

解决远程依赖问题

我们在上一节中看到，由于消失和爆炸梯度，香草 RNN 难以有效地学习远程依赖性。为了解决这个问题，通过 SEPP Hochreiter 和 JürgenSchmidhuber 于 1997 年开发了长短的短期内存网络。在 2014 年引入了门控经常性单位，并提供了更简单的 LSTM 版本。让我们回顾 LSTM 和 GRU 如何解决学习远程依赖性的问题。

长期记忆

LSTM 在每个时间步中引入了其他计算。但是，对于时间步长[ h [t] ），返回状态为内部状态[ f [t] ），并将其转发到下一个时间步。但是，在内部，这些向量的计算方式不同。 LSTM 引入了三个新的门：输入（o[t]），忘记（g[t]）和输出（c[t]）门。每个时间步也具有内部隐藏状态（ h [ t ] ）和内部存储器（σ(W x[t] + U s[t-1])。这些新单位的计算方式如下：

现在，让我们了解一下计算。门o[t]，σ(·)和[0, 1]是通过将其值限制在f[t]内的 S 型激活h[t]生成的。因此，它们在与另一个变量相乘时仅散发值的一小部分而充当门。输入门o[t]控制要保留的新计算输入的分数。遗忘门c[t]确定前一个时间步的效果，输出门f[t]控制要释放的内部状态量。内部隐藏状态是根据输入到当前时间步以及上一个时间步的输出来计算的。请注意，这与在普通 RNN 中计算内部状态相同。 h[t]是当前时间步的内部存储单元，它考虑了上一步的存储器，但缩小了一部分s[t]，并减小了内部隐藏状态的影响，但与输入门o[t]混合在一起。最后，z[t]将传递到下一个时间步，并根据当前内部存储器和输出门r[t]进行计算。输入，忘记和输出门用于有选择地包括以前的内存和当前的隐藏状态，这些状态的计算方式与普通 RNN 中的相同。 LSTM 的门控机制允许从较长的时间步长进行内存传输。

门控循环单元

GRU 比 LSTM 简单，只有两个内部门，即更新门（ z [t] ）和复位门（ r [t] ）。更新和重置门的计算如下：

*z[t] = σ(W<sup class="calibre28">z</sup>x[t] + U<sup class="calibre28">z</sup>s[t-1])*

*r[t] = σ(W<sup class="calibre28">r</sup>x[t] + U<sup class="calibre28">r</sup>s[t-1])*

使用输入 x [t]，状态 s [t- 1 上一个时间步中的]，更新和复位门：

由 S 形函数计算的更新确定在当前时间步中要保留多少前一步的内存。复位门控制如何将先前的存储器与当前步骤的输入合并。

与具有三个门的 LSTM 相比，GRU 具有两个门。它没有输出门和内部存储器，这两个都存在于 LSTM 中。 GRU 中的更新门决定了如何将先前的内存与当前的内存结合，并结合 LSTM 的输入门和忘记门所实现的功能。结合了先前存储器和当前输入的效果的复位门直接应用于先前存储器。尽管在沿序列传输内存的方式上存在一些差异，但 LSTM 和 GRU 中的门控机制都旨在学习数据中的远程依赖性。

使用哪一个-LSTM 或 GRU？

LSTM 和 GRU 都能够处理长 RNN 上的内存。但是，一个常见的问题是使用哪个？长期以来，LSTM 一直是语言模型的首选，这从它们在语言翻译，文本生成和情感分类中的广泛使用可以明显看出。与 LSTM 相比，GRU 的显着优势是可训练的砝码更少。它已应用于以前 LSTM 占主导地位的任务。但是，实证研究表明，在所有任务中，两种方法都无法胜过另一种方法。诸如隐藏单元的维数之类的模型超参数的调整可改善两者的预测。一个普遍的经验法则是在训练数据较少的情况下使用 GRU，因为它需要较少的可训练砝码。在大型数据集（例如用于开发语言翻译模型的数据集）的情况下，LSTM 被证明是有效的。

递归神经网络用于时间序列预测

我们将继续使用空气污染数据集来演示递归神经网络以进行时间序列预测。 LSTM 用于预测气压，而 GRU 在pm2.5上进行了演示。

读取和预处理数据与在 MLP 上的示例相同。原始数据集分为两组训练和验证，分别用于模型训练和验证。

makeXy函数用于生成回归变量和目标数组-X_train，X_val，y_train和y_val。由makeXy函数生成的X_train和X_val是形状为(number of samples, number of timesteps)的 2D 数组。但是，RNN 层的输入必须为(number of samples, number of timesteps, number of features per timestep)形状。在这种情况下，我们只处理pm2.5，因此number of features per timestep是一个。 Number of timesteps为七个，number of samples与X_train和X_val中的number of samples相同，它们被重塑为 3D 阵列：

X_train, X_val = X_train.reshape((X_train.shape[0], X_train.shape[1], 1)), 
                 X_val.reshape((X_val.shape[0], X_val.shape[1], 1)) 
print('Shape of 3D arrays:', X_train.shape, X_val.shape)

X_train和X_val已重塑为 3D 阵列，并且它们的新形状在前面的print语句的输出中可见，如下所示：

Shape of 3D arrays: (35057, 7, 1) (8753, 7, 1)

要将 LSTM 层添加到神经网络，我们需要从 Keras 导入 LSTM 类：

from keras.layers.recurrent import LSTM

用于开发时间序列预测模型的神经网络具有输入层，该输入层馈入 LSTM 层。 LSTM 层有七个时间步，这与为第二天的气压做出预测所用的历史观测值相同。仅 LSTM 的最后一个时间步返回输出。 LSTM 层的每个时间步中都有 64 个隐藏的神经元。因此，LSTM 的输出具有六十四个功能：

input_layer = Input(shape=(7,1), dtype='float32') 
lstm_layer = LSTM(64, input_shape=(7,1),  return_sequences=False)(input_layer)

接下来，将 LSTM 的输出传递到一个 dropout 层，该 dropout 层在传递到 output 层之前会随机丢弃 20％的输入，后者具有一个具有线性激活函数的隐藏神经元：

dropout_layer = Dropout(0.2)(lstm_layer) 
output_layer = Dense(1, activation='linear')(dropout_layer)

最后，将所有层都包裹在keras.models.Model中，并使用 Adam 优化器训练 20 个纪元以最小化 MSE：

ts_model = Model(inputs=input_layer, outputs=output_layer) 
ts_model.compile(loss='mean_squared_error', optimizer='adam')

我们已将keras.callbacks.ModelCheckpoint用作回调，以在验证集上跟踪模型的 MSE，并保存来自提供最小验证错误的时期的权重。验证集的最佳模型的 R 平方是 0.9959。下图显示了前 50 个实际值和预测值：

图 5.20：实际和 LSTM 预测的气压时间序列

该示例的代码可以在 Jupyter 笔记本code/Chapter_5_Air Pressure_Time_Series_Forecasting_by_LSTM.ipynb中找到。

我们使用了两个堆叠的 GRU 层来开发基于 RNN 的pm2.5时间序列预测模型：

gru_layer1 = GRU(64, input_shape=(7,1), return_sequences=True)(input_layer) 
gru_layer2 = GRU(32, input_shape=(7,64), return_sequences=False)(gru_layer1)

第一个 GRU 从先前的输入层获取顺序输入。第一个 GRU 的每个时间步都返回一个 64 维特征向量作为输出。该序列作为输入传递到下一个 GRU 层。第二个 GRU 层仅从最后一个时间步返回输出。使用 Adam 优化器对神经网络进行训练，以最大程度地降低 MAE 损失。

验证集获得的最佳 MAE 为 11.388。

该示例的完整代码在 Jupyter 笔记本code/Chapter_5_PM2.5_Time_Series_Forecasting_by_GRU.ipynb中。

图 5.21：实际和 GRU 预测的气压时间序列

卷积神经网络

本节描述了卷积神经网络（CNN），当输入数据为图像时，它们主要用于开发有监督和无监督的模型。通常，将二维（2D）卷积应用于图像，但是可以将一维（1D）卷积应用于图像上。顺序输入以捕获时间依赖性。本节将探讨这种方法，以开发时间序列预测模型。

2D 卷积

让我们从描述 2D CNN 开始，我们将导出 1D CNN 作为特殊情况。 CNN 利用图像的 2D 结构。图像的矩形尺寸为w，其中n为高度，hxwxn为图片的宽度。每个像素的颜色值将成为模型的输入特征。使用具有 28 x 28 个神经元的完全连接的密集层，可训练的权重数将为 28 x 28 x 100 =78400。对于 MNIST 数据集中的手写数字 32 x 32 的图像，第一个密集的可训练权重数具有 100 个神经元的层将为 c =3。CIFAR-10 数据集通常用于训练对象识别模型。此数据集中的彩色图像为 32 x 32 x 3 x 100 = 307200，并具有mxmx33 个颜色通道，红色，绿色和蓝色。因此，具有 100 个神经元的完全连接的密集层将具有可训练的权重。因此，训练图像的密集层成为计算难题。

CNN 通过将神经元仅连接到图像的局部斑块来解决此问题。如下图所示，尺寸过滤器应用于局部图像补丁。滤镜的第三维与图像的颜色通道数相同。滤波器中每个神经元的权重乘以图像中相应的像素值。该局部补丁的最终特征是通过将这些单个值相加，并可选地加上一个偏差，然后将总和传递给激活函数来计算的。激活功能的常用选择是整流线性单位（ReLu）：

ReLu（z）= 0 且 f z≤0

如果 z > 0 ，则 = z

ReLu 的一阶导数为 0 或 1，并且具有良好的梯度流动特性。因此，在训练具有多层的深度卷积网络中是优选的。下图所示的滤波器是一个卷积层。为了学习不同的特征，卷积层通常具有多个过滤器。

图 5.22：滤波器在卷积神经层生成特征中的应用

为了覆盖整个图像，以 4 x 4 像素单位的水平步幅和 2 x 2 像素单位的垂直步幅移动滤镜。滤镜在图像上的运动称为卷积。由整个图像上的卷积产生的特征形成一个矩形特征图。现在，让我们了解卷积如何生成特征图。为此，让我们考虑具有单色通道并且像素值如下图所示的图像。我们应用 3 x 3 卷积滤波器，并将其水平移动一个像素单位，垂直移动一个像素单位。

第一卷积被计算为 1×-1 + 2×1 + 2×2 + 1×1 = 6。其他卷积的特征也以类似的方式计算。在此示例中，所有计算值均为正，ReLu 激活的作用就像一个身份函数。这个卷积过程将创建一个特征图，如下图所示：

图 2.23：二维卷积图

请注意，2 x 2 滤镜的单位像素水平和垂直跨度将特征图的大小从原始图像中缩小为3x3。6x6。为了确保特征图的大小与图像的大小相同，我们可以在图像的水平和垂直边界上串联零值像素，并使用滤镜。添加零值像素称为零填充，在这种情况下，零填充将使输入图像的尺寸为 2 x 2。下图显示了从零填充图像生成的特征图：

图 2.24：二维卷积生成的特征图

由卷积层生成的特征图将被馈送到下游的卷积层，就像原始输入已被提供给第一卷积层一样。彼此叠置的多个卷积层从原始图像生成更好的特征。然后将这些特征传递到下游的全连接密集层，这些密集层在对象类集合上生成 softmax 输出。

在大多数基于图像的深度学习模型中，不使用卷积层对原始图像或中间特征图进行下采样。如果是这种情况，将卷积层的输出馈送到密集层时，可训练权重的数量仍然很大。那么卷积层能实现什么呢？我们如何对特征图进行下采样以减少密集层中可训练权重的数量？

首先，尽管我们可以在卷积层中实现下采样，但最好将它们用于提取图像特征，例如边缘，拐角，形状等。其次，由池化层实现的下采样，该池化层将过滤器应用于要素图的本地面片并计算单个要素。过滤器卷积在整个特征图上。池化层的卷积生成一个降采样的特征图。池化层没有可训练的权重，但是应用简单的算术函数（例如最大值或平均值）来生成其输出特征图。实证研究表明，maxpooling 层可以提取有用的图像特征并提高图像识别模型的准确性。下图说明了 2 x 2 maxpooling 层如何生成特征图。 maxpooling 层的水平和垂直步幅为两个单位，并且已将输入特征图从降低采样到 28 x 28。

图 2.25：二维 maxpooling 的图示

下图说明了用于图像识别的神经网络。有两个卷积层块，然后是一个完全连接的密集层块。带有标签的图像被馈送到该网络，并且输出是标签类别上的概率。例如，可以使用具有t个手写数字（0、1、2，... 9）灰度图像的 MNIST 数据集训练数字识别模型。灰度图像被馈送到神经网络，如下图所示。网络的最后一层是 softmax 层，它提供了十位数字类的预测概率。概率最高的类别是预测数字。

图 2.26：用于图像分类的深度卷积神经网络架构

一维卷积

一维卷积层可用于开发时间序列预测模型。具有 1 个m观测值的时间序列类似于尺寸为p的图像，其高度为单个像素。在这种情况下，可以使用1x3滤波器将 1D 卷积作为 2D 卷积的特殊情况。另外，过滤器仅沿水平方向移动1x8时间单位的步长。

让我们了解一维卷积的工作原理。考虑下图，该图显示了十个时间步长的时间序列。（1 x -1）+（2x1）+（-1 x 2）= -1 过滤器在序列上移动了一个时间单位。因此，生成了1x3特征图。特征图的第一元素被计算为1x10。其余时间步长以类似的方式计算。我们没有对原始时间序列进行零填充，因此特征图比原始时间序列短两个单位。但是，对输入时间序列的开始和结束进行零填充并使用相同的滤波器将导致1x3特征图。可以将池化层与卷积层堆叠在一起以对特征图进行下采样。

图 2.27：一维卷积图

使用 1 x 3 卷积滤波器的方法等效于训练几个三阶的局部自回归模型。这些局部模型在输入时间序列的短期子集中生成特征。在 1D 卷积层之后使用平均池化层时，它将在前一卷积层生成的特征图上创建移动平均值。此外，几个 1D 卷积和合并层相互堆叠时，提供了一种从原始时间序列中提取特征的有效方法。因此，在处理复杂的非线性时间序列（例如音频波，语音等）时，使用 CNN 证明是有效的。实际上，CNN 已成功应用于音频波的分类。

在以下部分中，我们将使用 1D 卷积开发时间序列预测模型。

一维卷积用于时间序列预测

我们将继续使用空气污染数据集来演示时间序列预测的 1D 卷积。输入到卷积层的形状为(number of samples, number of timesteps, number of features per timestep)。在这种情况下，number of timesteps是七个，number of feature per timestep是我们仅担心的空气压力，这是一个单变量。要开发卷积神经网络，我们需要导入三个新课程：

from keras.layers.convolutional import ZeroPadding1D 
from keras.layers.convolutional import Conv1D 
from keras.layers.pooling import AveragePooling1D

在输入层之后添加ZeroPadding1D层，以在每个序列的开头和结尾添加零。零填充可确保卷积层不会减小输出序列的维数。在卷积层之后添加的池化层用于对输入进行下采样：

zeropadding_layer = ZeroPadding1D(padding=1)(input_layer)

接下来，我们添加一个Conv1D层。 Conv1D的第一个参数是过滤器的数量，它们确定输出中要素的数量。第二个参数指示一维卷积窗口的长度。第三个参数是strides，表示移动卷积窗口的位数。最后，将use_bias设置为True会在计算输出特征时增加一个偏差值。在这里，一维卷积可以认为是在三个时间单位的滚动窗口上生成本地 AR 模型。

conv1D_layer = Conv1D(64, 3, strides=1, use_bias=True)(zeropadding_layer)

接下来添加AveragePooling1D，通过对三个时间步进行平均（跨度为一个时间步）对输入进行下采样。在这种情况下，平均池可被视为对三个时间单位的滚动窗口进行移动平均。我们使用平均池而不是最大池来生成移动平均值：

avgpooling_layer = AveragePooling1D(pool_size=3, strides=1)(conv1D_layer)

前面的池化层返回 3D 输出。因此，在传递到输出层之前，先添加Flatten层。 Flatten层将对(number of samples, number of timesteps×number of features per timestep)的输入进行整形，然后在通过Dropout层之后将其馈送到输出层：

flatten_layer = Flatten()(avgpooling_layer) 
dropout_layer = Dropout(0.2)(flatten_layer) 
output_layer = Dense(1, activation='linear')(dropout_layer)

将这些层打包到keras.models.Model包装器中，并使用 Adam 优化器对其进行了培训，以使 MAE 降至最低。最佳模型的验证 R 平方为 0.9933。

详细的实现在 Jupyter 笔记本code/ Chapter_5_Air_Pressure_Time_Series_Forecasting_by_1D_Convolution.ipynb中。我们采用了类似的方法来为pm2.5开发基于 1D CNN 的时间序列预测模型，并将代码编写在 Jupyter 笔记本code/ Chapter_5_PM2.5_Time_Series_Forecasting_by_1D_Convolution.ipynb中。

概括

在本章中，我们描述了三种基于深度学习的方法来开发时间序列预测模型。神经网络适用于几乎没有诸如长期趋势和季节性之类的基本属性的信息，或者这些信息过于复杂而无法通过传统统计方法以可接受的准确度进行建模的情况。诸如 MLP，RNN 和 CNN 之类的不同神经网络架构会从数据中提取复杂的模式。如果使用适当的措施对神经网络模型进行了训练，以避免对训练数据过度拟合，那么这些模型可以很好地针对看不见的验证或测试数据进行概括。为了避免过度拟合，我们应用了 dropout，它在深度神经网络中广泛用于各种数据集和应用。我们希望本章能使您对可用于时间序列预测的高级技术有所了解。本章随附的 Jupyter 笔记本有望为您提供必要的基础知识，这对于开发用于解决复杂问题的更高级模型很有用。

六、Python 入门

当您选择阅读本书时，我们认为您可能已经掌握 Python 的使用知识，如果不是，那是一位亲自动手的专家，他恰巧可以生存和呼吸 Python。如果至少您对 Python 有一定的了解，则可以选择跳过此附录。如果您不熟悉 Python 或正在寻找如何开始使用编程语言，阅读本附录将帮助您克服最初的障碍。这也将为您带来享受本书各章所需的知识。因此，事不宜迟，让我们开始吧！

Python 是一种通用的高级解释型编程语言，于 1991 年问世。它的创建者 Guido van Rossum 于 1989 年圣诞节开始编写该语言的解释，并以他最喜欢的电视之一命名该语言。展示了 Monty Python 的飞行马戏团。

Python 强调通过空格缩进来界定代码块而不是大括号，从而强调代码的可读性，而在 C，C ++ 和 Java 中，大括号被广泛使用。 Python 的另一个强大功能是它的简洁性，它使程序员可以用比 C，C ++ 和 Java 更少的代码行来表达概念。例如，在 Python 开发人员中，使用 lambda 函数作为仅一行声明函数的一种快速有效的方法就很受 Python 开发人员的欢迎。（不用担心；稍后我们将介绍 lambda 函数。）除这些以外，Python 支持动态类型绑定并自动为您处理内存。此外，Python 解释器可用于所有主要操作系统。 CPython 是解释器的最流行的开源实现。此外，Python 受到数百种广泛功能包的支持，这些功能包包括网络开发，GUI 构建，高级内存管理，用于多种文件格式的文件处理，科学和数值计算，图像处理，机器学习，深度学习，大数据，和许多其他。 PyPI 是 Python 软件包的官方存储库，几乎所有知名软件包都可以在那儿找到。指向 PyPI 网站的链接为这个页面。但是，在安装之前，我们不必从此网站下载软件包。有命令行（CL）工具可以为我们完成工作。一旦安装了该语言的基本版本，即可使用这些工具。而且，通过确保所请求的软件包与当前版本的 Python 兼容并且已经安装或需要安装依赖项，这些 CL 工具需要我们大量工作。

本附录涵盖以下主题：

安装
基本数据类型
关键字，控制语句和功能
迭代器和生成器
类和对象

安装

在您的计算机上设置 Python 需要安装三个组件-解释程序，集成开发环境（IDE）和支持应用开发的软件包。值得庆幸的是，我们不必单独安装它们。有将这三个捆绑在一起的软件包，并且一次安装可以最大程度地利用您的计算机。

Python 安装程序

有两种广泛使用的安装程序选项，它们将解释器，IDE 和有用的 Python 软件包捆绑在一起：

选项 1：捆绑了所有三个组件的安装程序可从这个页面下载。该网站提供了几种不同版本的语言的安装程序。在编写本书时，最新版本是 Python 3.2.6 和 Python 2.7.13。安装程序附带一个基本的 IDE，是编写一些代码并了解该语言工作方式的好地方。这些安装程序提供了一些软件包，例如文件处理，内存管理，数学计算等。但是，有关科学计算，统计建模和机器学习的大多数软件包都需要与它们的依赖项一起单独安装。这使此选项很耗时。
选项 2：我们可以从这个页面下载免费的安装程序。 Continuum Analytics 在将解释器，IDE 和数百个有用的软件包组合到同一安装程序中的工作非常出色。这些软件包涵盖了您在数据科学项目上工作可能需要的大部分内容。捆绑包中还包括 Jupyter Notebook 和 Jupyter Lab，这是流行的基于 Web 的 Python 开发工具。软件包的完整列表位于这个页面。这些安装程序被称为 Python 的 Anaconda 发行版。 Python 2.7 和 3.6 版本均可用。本书示例中使用的大多数软件包都是 Anaconda 发行版的一部分。我们只需要分别安装 Keras 和 Tensorflow。

下一节将介绍本书示例的所需步骤。

运行示例

本书的代码示例是使用 Jupyter Notebook 开发环境编写的。所有代码文件都具有.ipynb 扩展名，可以在该书的 GitHub 存储库的代码文件夹中找到。

要运行 Jupyter Notebook，您需要安装 Anaconda Python Distribution。我们从这个页面下载了 Python 3.6 版本的发行版。

下载安装程序后，需要执行以下步骤来为示例设置环境：

安装 Anaconda Python 发行版。按照标准安装的提示进行操作。安装程序会将 Python 可执行文件的路径添加到 PATH 环境变量中。安装程序会提示您添加到 PATH 变量。阅读说明后，请确保单击“是”或“确定”，这将允许安装程序更新 PATH。
打开一个新的命令提示符窗口，然后运行以下命令：

python --version

这是为了确保将 Python 可执行文件添加到 PATH 变量中。
现在是时候设置编程环境来运行 Jupyter Notebook 了：
- 将图书的 GitHub 存储库“实用时间序列”克隆或下载到本地磁盘，例如 D。
- 打开一个新的命令提示符，然后运行以下命令：

D: (optional if default location is some other drive)
cd "Practical Time Series"
jupyter notebook -notebook-dir="."

第三步启动 Jupyter Notebook 服务器，并在默认浏览器中自动打开主页。主页显示磁盘文件夹D:\Practical Time Series的文件夹和文件。单击主页中的code文件夹，您将看到 Jupyter 笔记本。每个笔记本都有以下行，用于设置笔记本的当前工作目录：

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

如果已将Practical Time Series文件夹保存在 D 驱动器中，则无需更改此行。否则，必须更改它以反映实际的当前工作目录。例如，如果共享文件夹另存为C:\Local\Practical Time Series，则将前一行代码更改为以下内容：

os.chdir('C/Local/Practical Time Series')

Jupyter 笔记本中的代码段显示在带有灰色背景的矩形单元格中，如以下屏幕截图所示。要运行笔记本，请单击“运行单元”按钮（按钮面板左侧的第八个按钮），从顶部开始依次执行每个单元：

基本数据类型

像所有主要的编程语言一样，Python 支持基本的数字数据类型，例如 int，long 和 float。 None 表示 Python 中的空指针。

字符串是 Python 中的顺序数据类型。 Python 中常用的其他顺序数据类型是列表和元组。列表和元组之间的区别在于，前者是可变的，而后者是不可变的类型。因此，如果您尝试修改元组，则解释器将引发错误。让我们更深入地研究列表和元组。

列表，元组和集合

列表是元素的集合。元素可以是任何数据类型，列表可以包含不同数据类型的元素。列表具有几个重要功能，例如append，extend，insert，index，pop等。下表总结了这些功能的功能：

列表还具有其他一些功能。有关更全面的摘要，我们建议使用任何教授 Python 语言的材料。

上表中提到的功能在 Jupyter Notebook code/Getting_started.ipynb中进行了演示。请注意，Python 使用零作为序列第一个元素的索引。因此，mylist[0]是第一个元素，mylist[1]是第二个元素，依此类推。

与列表不同，元组是不可变的序列，并且在声明后就不允许进行任何更改。元组仅支持两个函数，即 count 和 index，它们返回元组中元素的数量和现有元素的索引。请求不存在的元素的索引将引发ValueError。

通过将变量作为参数传递给函数，len函数可用于查找列表的长度以及元组。

此外，列表和元组是可迭代的。可迭代的技术定义是返回迭代器的任何对象。在以下各节中，我们将详细讨论可迭代器和迭代器。但是，暂时请记住，迭代器是对象，可以通过为对象运行循环来评估其各个元素。以下代码演示了对列表和元组的迭代：

#iterating over a list
print('Elements in mylist:')
for i in mylist:
print(i)
print('Elements in mytuple:')
for i in mytuple:
print(i)

前面的代码段的输出如下：

Elements in mylist:
-1
2
3
4
5
6
[7, 8]
7
8
Elements in mytuple:
1
2
3
4

列表和元组中的元素也可以切成一定范围的索引。切片是从原始序列创建子集的有效方法：

#slicing lists and tuple
print('Second to fifth elemnts of mylist:', mylist[1:6])
print('Last three elemnts of mytuple:', mytuple[-3:])

这将返回以下输出：

Second to fifth elements of mylist: [2, 3, 4, 5, 6]
Last three elements of mytuple: (2, 3, 4)

请注意，在前面的代码中用于切片mytuple的负索引。负索引允许从末尾开始遍历。因此，-3指示从末端开始的第三个位置。从末尾开始的第一个位置是-1的索引。

集是唯一元素的有序序列，并支持标准集操作，例如联合，交集和差异。通过传递可迭代对象来创建集合。集合支持弹出操作，但与列表不同，弹出功能仅删除集合的第一个成员。在 Jupyter 笔记本code/Getting_started.ipynb中给出了集合的示例。

弦乐

字符串在 Python 2.x 中是 8 位字符的可迭代，在 Python 3.x 中是 Unicode 字符。可以用多种方式声明字符串，以下代码片段中显示了其中的一些方式：

#string declarations
mystring = 'Practical Time Series'
print(mystring)
mystring = "Practical Time Series"
print(mystring)
mystring = """Practical Time Series"""
print(mystring)
mystring = '''Practical Time Series'''
print(mystring)

前面的所有四个声明均导致相同的字符串，并且print函数在每种情况下均提供相同的输出：

Practical Time Series

字符串是 class，str 类型的对象，并支持使字符串处理变得容易的几个函数。我们在 Jupyter 笔记本code/Getting_started.ipynb中显示了一些示例。

地图

Python 支持哈希图或关联数组。哈希图中的每个元素都与定义或键相关联。 Dict 是 Python 中的基本 hashmap 实现。字典中的元素称为值，它们映射到称为键的引用。可以从键值对的可迭代实例化 dict。这个可迭代的对象可以是列表，元组或迭代器，我们将在接下来的部分中进行讨论。以下代码段显示了创建字典的多种方法：

mydict = dict([(0,1), (1,2), (2,3), (3,4), (4,5)])
mydict = dict([[0,1], [1,2], [2,3], [3,4], [4,5]])
mydict = {0:1, 1:2, 2:3, 3:4, 4:5}

可以通过调用将键作为输入的 get 方法来检索字典中的元素。另一种方法是使用dictionary_variable[key]。对于不存在的键，get 函数返回 None，但是第二种方法抛出KeyError错误。函数项返回键-值对的迭代器。所有对都不会在列表中返回，而是可以在迭代器上循环。

Python 是一种面向对象的编程语言，因此支持诸如对象的复杂数据类型，这些对象是类的实例，并具有属性和方法。以下各节之一将介绍类和对象。

关键字和功能

关键字是保留字，不能用作变量名。下表列出了关键字及其用途：

在 Python 中，函数分为三种类型：

内联函数
正常功能
内建功能

内联函数使用 lambda 关键字声明。这种类型的函数声明可以命名或匿名。匿名内联函数用于在另一个函数或代码块内声明一个函数，该函数或代码块循环一个序列并将内联函数应用于其元素。命名的内联函数是单行函数声明，但具有函数名称：

#Inline anonymous function
mylist = [0, 1, 2, 3, 4, 5]
processed_list = list(map(lambda x: x+1, mylist))

#Inline named function
myfunc = lambda x: x+1
processed_list = [myfunc(i) for i in mylist]
print(processed_list)

请注意，在前面的代码块中使用了map函数。 map是 Python 中的预定义函数之一。它接受一个匿名内联函数和一个 Iterable 作为输入，并将该内联函数应用于 Iterables 的每个元素。 map函数的输出是一个新的可迭代函数，具有针对输入序列的每个元素的，由内联函数返回但由内联函数转换的值。

普通功能是使用 def 关键字定义的。如果没有return语句，普通函数将在其范围内执行代码并返回 None。使用 return 语句，可以在调用函数时在函数外部获取普通函数的输出：

#Normal function which returns None
def myfunc(a):
print(Within the function:', a+1)
print('Return value of the function:', myfunc(5))

打印功能的输出如下：

Within the function: 6
Return value of the function: None

#Normal function which returns an output
def myfunc(a):
a = a+1
print('Wihtin the function:', a)
return a
print('Return value of the function:', myfunc(5))

最后一个打印功能的输出如下：

Within the function: 6
Return value of the function: 6

最后，Python 具有几个预定义的功能，例如range，map，filter，zip，enumerate等。富有表现力和简洁是 Python 设计哲学的两个关键要素。为此，预定义功能仅需几行代码就可以帮助我们实现更多目标。下表描述了预定义功能的功能：

在 Python 2.x 中，所有这些函数都将列表作为输出返回，但是在 Python 3.x 中，它们返回迭代器。为了获得实际的输出，必须在迭代器上循环。如果输入序列太大而又太大而无法容纳计算机的内存，则此设计很有用。例如，考虑一个巨大的文本文件，它必须与行号一起逐行读取。在 Python 中，可以通过遍历文件读取器返回的迭代器逐行读取文件。为了获得行号以及行号，我们可以使用文件迭代器调用枚举。枚举函数将文件迭代器视为一个序列，并生成一个（line_number，行）对的迭代器，可以将其逐个循环，而无需将整个文件加载到内存中。我们刚刚提到的其他内置函数也可以采用类似的方法。在 Jupyter Notebook code/Getting_started.ipynb中，我们展示了如何将这些功能与文件阅读器一起使用。我们建议您对其他类型的序列尝试类似的方法。

我们已经介绍了一些最常用的内置函数。但是，建议您参考有关 Python 编程语言的书或教程以获取详尽的列表。

迭代器，可迭代对象和生成器

在 Python 中，我们经常遇到迭代器，可迭代器和生成器，因为它们是循环作为序列的数据类型或数据结构或可以从中创建序列的有效方法。使用这些循环技术的一个明显优势是它们需要更少的内存。因此，当您必须逐个元素地访问序列时，这些技术将变得非常有用，因为不需要将大序列立即全部加载到内存中。例如，如果您需要找到前一万亿个正整数的平方，则无需创建一个数据结构来同时将所有数字保存在内存中。迭代器，可迭代器和生成器可用于顺序生成和处理这些数字。另一个示例是处理大型文本文件。整个文件可能不适合内存。因此，如果我们需要处理文件，例如查找文件每行的字数，则可以迭代遍历这些行并逐行处理它们。由于迭代器，可迭代对象和生成器是如此有用，因此让我们了解它们是什么以及如何使用它们。

迭代器

迭代器是具有__iter__和 next 函数的类的实例的对象，并且可以与 for 循环一起使用以逐个顺序地遍历序列。 __iter__函数可将对象识别为迭代器。调用下一个函数以逐个获取序列的元素。每次调用下一个函数时，它都会从预定义的序列中返回一个元素，或者创建一个元素。因此，下一个功能可以实现创建元素的逻辑。没有更多元素时，__next__函数将引发 StopIteration 错误。

让我们通过示例来了解迭代器是如何工作的。我们创建一个迭代器以从预定义的序列中返回元素：

class MyIterator(object):

def __init__(self, seq):
self.seq = seq
self.i = 0

def __iter__(self):
return self

def next(self):
if self.i < len(self.seq):
i = self.i
self.i += 1
return self.seq[i]
else:
raise StopIteration()

创建 MyIterator 类的对象，然后调用下一个函数五次。前四个调用从序列中返回一个整数，但是最后一个调用导致下一个函数抛出StopIteration错误，因为到目前为止已经遍历了序列的整个长度：

itr = MyIterator([1,2,3,4])
print(itr.next())
print(itr.next())
print(itr.next())
print(itr.next())
print(itr.next())

现在让我们声明一个迭代器，该迭代器在下一个函数中实现数据生成逻辑，而不是从预定义序列中返回元素：

class MyIterator(object):

def __init__(self, n):
self.n = n
self.count = 0

def __iter__(self):
return self

def next(self):
if self.count < self.n:
i = self.count
self.count += 1
return i
else:
raise StopIteration()

创建MyIterator类的对象，然后调用下一个函数六次。前五个调用返回下一个函数生成的正整数，但最后一次调用引发 StopIteration 错误：

itr = MyIterator(5)
print(itr.next())
print(itr.next())
print(itr.next())
print(itr.next())
print(itr.next())
print(itr.next())

可迭代

不具有__iter__或 next 函数但可用于创建迭代器的对象是可迭代的。内置的__iter__函数将顺序对象作为输入并返回迭代器。然后，另一个内置函数接下来将使用迭代器，并在每次调用时返回迭代器的元素。列表和元组不是迭代器，但可用于使用 next 函数创建迭代器。以下代码片段演示了如何使用列表作为迭代器：

mylist = [1,2,3,4]
mylist_iter = iter(mylist)
print(type(mylist))
print(type(mylist_iter))

<class 'list'>
<class 'list_iterator'>

请注意，mylist的类型为 list，而由 iter 函数创建的mylist_iter的类型为list_iterator。

发电机

生成器提供了迭代器的功能，但无需编写整个类。在许多情况下，创建序列元素的逻辑是在迭代器的下一个功能中实现的。例如，需要在迭代器的下一个函数中编写用于逐行读取文本文件并返回每一行中特定单词的出现的代码。仅在满足开发迭代器的要求的情况下，在该类中实现了其他功能，例如__init__和__iter__。

生成器是简化开发的特殊功能。函数中的 yield 关键字表示 Python 解释器该函数是生成器。使用生成器函数，我们仅实现在 while 循环中在生成器函数中创建序列元素的逻辑。每次调用函数时，yield 函数都会返回元素。

让我们实现一个生成器函数，该函数在每次被调用时返回一个整数：

def int_gen():
count = 0
while True:
i = count
count += 1
yield i

我们将生成器函数分配给一个变量，并在该变量上调用内置的 next 函数。下一个函数将变量作为参数运行五次：

ig = int_gen()
print(next(ig))
print(next(ig))
print(next(ig))
print(next(ig))
print(next(ig))

打印输出如下：

每次在 ig 上调用下一个函数时，它都会返回 count 变量的当前值并将其递增 1。计数变量保持在生成器的内部状态，因此可以返回其最新值。

类和对象

类是变量和函数的逻辑分组。 class关键字在 Python 中用于定义此类逻辑分组。类通常代表现实生活中的实体，例如，书籍，作者，出版商等。实体具有属性，这些属性由类中定义的变量表示。类中的函数（通常称为方法）定义了如何捕获和转换有关实体实例的数据。类的实例是实体的单个实现。例如，书是一个实体，而实践时间序列分析是书的一个实例。要创建实例，我们初始化一个类的对象。对象定义涉及通过构造函数将值分配给类的变量。该工作由__init__方法完成，该方法接受输入并将其分配给类变量。 __init__方法可以基于创建对象的逻辑在内部调用其他函数。让我们定义一个关于书籍的类：

import datetime
class Book(object):

def __init__(self, name, date_of_publication, nb_pages, publisher):
self.name = name
self.date_of_publication = datetime.datetime.strptime(date_of_publication, '%Y-%m-%d')
self.nb_pages = nb_pages
self.publisher = publisher
self.authors = []

def add_author(self, author_name):
self.authors.append(author_name)

def print_date_of_publication(self, print_format='%Y-%m-%d'):
print(self.date_of_publication.strftime('%d-%m-%Y'))

现在，我们将创建Book类的对象：

mybook = Book('Practical Time Series Analysis', '2017-09-15', 200, 'Packt')

请注意，__init__方法未指定书的作者，而是创建了一个空列表，即作者。我们可以在mybook对象上调用 add_author 函数以添加作者。另外，date_of_publication最初设置为%Y-%m-%d格式。 print_date_of_publication函数以 print_format 作为输入，并以另一种格式显示日期。例如，我们可以将发布日期打印为%d-%m-%Y：

mybook.print_date_of_publication(print_format='%d-%m-%Y')

概括

本附录介绍了 Python 编程语言的基础。已经讨论了诸如数据类型，关键字，函数，类，迭代器，可迭代和生成器之类的主题。这些编程技术构成了 Python 的构建块。

本书的各章使用了此处已讨论的几种概念和编程技术。

posted @ 2025-10-27 08:55 绝不原创的飞龙阅读(12) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

掠夺·扩张·投机·博弈