深入研究空气质量数据
这里有一个按时间戳索引的数据集。您的数据可能涉及储备需求和供应，您的任务是预测战略产品的理想补充间隔。或者，您可能需要将历史销售信息转化为团队的关键可操作见解。也许您的数据是关于历史利率和一些股票价格的金融数据。也许您的任务是建模市场波动性，并需要量化投资周期内的货币风险。无论是从社会科学、能源分配到医疗保健再到环境研究，案例很多。但这些情景有什么共同点呢？第一，您手头有一个时间序列任务。第二，您肯定会从开始一个简洁而全面的探索性分析中受益。
目录
本文的目标
但是，执行探索性时间序列分析意味着什么？与其他数据科学问题不同，从时间序列数据中获取见解可能会棘手且远非一帆风顺。您的数据可能具有重要的基础趋势和季节性，或者适合在其复杂的周期模式中进行嵌套预测。区分由数据生成过程中的故障引起的异常离群值与持有关键信息的实际异常可能具有挑战性。而处理缺失值可能并不像您期望的那样简单。
本文将概述我在研究时间序列数据集时取得的成功过程。我将沿着我探索细颗粒物（也称为PM 2.5）的测量数据前进，这是空气污染和空气质量指数的主要成分之一。我将重点介绍一些最佳实践，特别关注详细生成清晰且高度信息化的可视化和统计摘要。
数据集描述
这里研究的数据来自加拿大不列颠哥伦比亚省温哥华市的四个监测站。它们记录了细颗粒物PM 2.5（直径小于或等于2.5微米的细小颗粒）的每小时平均测量结果，单位为µg/m3（每立方米的微克），日期从2016年1月1日到2022年7月3日。
PM 2.5主要来自化石燃料的燃烧，在城市中，它通常来自汽车交通和建筑工地。污染物的另一个主要来源是森林和草地火灾，它们很容易被风带走[1]。
下图显示了我们将要探索的站点的大致位置。
库和依赖项
我们将使用Python 3.9和绘图库Matplotlib以及Seaborn进行可视化。对于统计测试和数据探索，我们将使用statsmodels Python模块和SciPy库。所有的数据处理和辅助任务将使用Pandas和Numpy来处理。
这些软件包在流行的Python发行版和托管笔记本，如Anaconda和Miniconda、Google Collab或Kaggle Notebooks中都可以使用。因此，这里的每个代码示例在您选择的环境中都应该很容易复现。
入门指南
首先导入我们的库，我们将调用matplotlib.dates.mdates和datetime模块来帮助我们处理DateTime索引。为了生成一致的可视化效果，我还喜欢从定义我们的绘图样式和颜色调色板开始。接下来让我们开始。
这里是我们将深入研究的一个广泛视图——我们将花一些时间来对我们的数据有第一次了解。
为了进行这个可视化，我们将使用Seaborn的关系图，它将从我们DataFrame的聚合长版本中读取数据。为此，我们将使用Pandas的melt和resample方法，使用平均聚合在24小时间隔内。这将把我们的数据粒度从每小时减少到每日平均测量值，目的是减少生成图表所需的时间。
通过整个时间序列中所有四个监测站的清晰图片，我们已经可以开始做一些笔记：
- 存在一些主要异常情况，它们似乎在夏季和早秋时更为普遍。
- 这些异常情况似乎是由规模较大的事件引起的，因为它们在大致相同的时间段内影响了所有四个监测站。
- 如果你仔细观察，我在每个图表中都包含了所有站点的淡灰色散点图。通过这个细节，我们可以看到例如2017年的异常情况对北温哥华的两个站点产生了重大影响（它们达到了更高的PM 2.5值），而2018年的情况则相反。这种技术还确保了所有四个折线图在相同的Y轴范围内。
从这个初始图表中可以得出一些好的做法：
- Matplotlib允许高度可定制的轴刻度定位器和格式化器。在这个例子中，我使用了
mdates模块的MonthLocator创建了辅助月份定位器，有助于提高整体可读性。 - 确切的预期时间范围从我们的图表传递给了图表标题或副标题。“预期”是因为我们的可视化可能由于绘制期间两端的缺失值而被截断，这种做法可以帮助识别这些问题。这也是一种良好的文档实践，用于报告你的发现。
这是一个很好的开始，但让我们稍微缩小一下我们的视野。在接下来的章节中，我们将以更短的时间段来观察，但是现在使用我们最初的每小时数据。
详细观察
从现在开始，我们将开始定义一些可以快速调用以生成定制化可视化的函数。你可以把它看作是设置一个非常有帮助的分析工具集的手段，这将对我们前进的过程非常有帮助。
这个第一个函数将帮助我们查看特定时间段内的个体时间序列。我们将首先查看北温哥华Mahon Park站2017年的数据。
我们可以看到这里除了主要异常情况外，还存在一些局部较小的峰值。在2017年的年初和12月也出现了较高的波动性（方差在较短的时间跨度内增加）。
让我们进一步细看，超出异常范围，以便我们可以在Y轴的更窄范围内分析我们的数值。
我们可以看到这里有一些缺失值。我们函数中的fill=True参数可以帮助我们识别它们，并且这是一种突出显示数据缺失情况的好方法。那些最初很难注意到的小中断现在变得清晰可见。
你可能还注意到了X轴的自定义日期格式。对于上面的图表，我通过增加了定制的主要和次要定位器和格式化器来改进我们的plot_sequence()函数。这个新的功能可以根据可视化的时间跨度自适应我们的图表，并相应地格式化X轴。下面是包含在函数中的代码段：
现在我们知道我们的数据集有一些中断，所以让我们更仔细地研究一下。
缺失值
对于表格数据问题，在这一部分中，我们可能会着重定义MAR（随机缺失）和MNAR（非随机缺失）。然而，根据我们数据的性质（感官时间测量），我们知道数据流中的中断可能并非故意。因此，在这些情况下，更重要的是区分孤立的缺失值和连续的缺失值，以及完全缺失的样本和缺失值。这里的可能性很广泛，所以我将在将来专门为那个方面撰写一篇文章。
现在，让我们从查看缺失值热图（missing values heatmap）开始。我们将再次为此定义一个函数。
热图非常好，因为它们可以帮助我们量化缺失以及在时间轴上定位它们。从这里，我们可以得出以下结论：
- 我们没有完全缺失的样本（即缺失值在某个时间段内同时出现）。这是符合预期的，因为监测站点发送的数据是相互独立的。
- 我们的时间线早期有很长的缺失值序列，随着时间的推移，数据的可用性似乎有所改善。
某些后续部分的统计分析对于缺失值会有问题。因此，我们将使用简单的技术来处理它们，例如：
- Pandas的
ffill()和bfill()方法。它们分别用于向前或向后延伸到最近可用值。 - Pandas的线性或样条插值
interpolate()方法。它使用相邻的观测值来绘制曲线来填充缺失间隔。
间歇性
根据我们数据的性质，我们不应该预期负值。正如开始时提到的，我在数据预处理时将它们视为缺失值。让我们调用我们的摘要统计数据来确认这一点。
我们可以看到每个站点的最小测量值都是零，这引出了下一个问题。我们的时间序列是否是间歇性的？
当您的数据有大量恰好为零的值时，就具有间歇性。这种行为提出了特定的挑战，并且在模型选择过程中必须考虑到这一点。所以我们的数据中有多少零值？
我们可以看到零值的数量微不足道，因此我们没有间歇性系列。这是一个简单但至关重要的检查，特别是如果您的目标是预测。对于某些模型来说，预测绝对零可能很困难，并且如果您想预测需求的话，这可能是一个问题。例如，如果客户预期没有产品，那么您就不希望计划将三个产品交付给他。
季节性
了解时间序列中的周期对于规划建模阶段至关重要。如果您决定过度聚合数据，您可能会丢失较小周期的关键信息，或者它可以帮助您确定在较小颗粒度上进行预测的可行性。
我们将使用一些箱线图来开始研究这个问题。但首先，我们将暂时删除前5%的数据，以便我们可以更好地查看数据。
在这个下一个功能中，我们将使用一系列的箱线图来研究我们数据中的周期性。我们还将将颜色调色板映射到中位数值，以便它可以作为另一个有趣的视觉线索。
这个第一个图返回了每小时的测量值。在这里我们可以看到：
- 从上午9点到下午2点，PM 2.5的值一直较高。
- 北温哥华以外的站点在晚上8点到10点也显示出高峰。
- 清晨的值从凌晨2点到5点最低。
现在看看每周季节性和一周内不同值之间的差异。
从这里我们看到：
- 周末PM 2.5的值较低。
- 周二的污染水平稍微较高。
最后，看看月度趋势。
在这里，我们可以观察到：
- 8月份的PM 2.5值始终较高，各年都一样。
- 南方的站点在6月和7月的PM 2.5值较低，而北温哥华的站点在1月的测量值较低。
最后，从这些图中可以得出更多好的实践经验：
- 不要轻率地使用颜色调色板，因为它们可能会导致你产生错误的解释。如果我们只是将
pallette=”mako”传递给我们的箱线图，它将被映射到我们的X轴，而不是我们感兴趣的变量。 - 网格图在低维数据中是强大的信息容器，可以通过Seaborn的
relplot()或Matplotlib的subplots()快速设置。 - 您可以使用Seaborn的
boxplot()的order参数来按有意义的顺序重新排序X轴。我用它来重新组织星期几的X标签。
从我们的时间序列中可以得到更多关于季节性的详细信息的观察是可能的。然而，这将留给将来一篇关于时间序列相似性和模型选择的文章进行更深入的探讨。如果这是您感兴趣的话题，请务必在VoAGI上关注我以接收我的未来文章。
现在，让我们快速查看一下我们众所周知的统计系数之一，以探究我们四个站点之间的线性关系。
皮尔逊相关系数
R程序员可能熟悉下面的图。一个相关矩阵图是一种简洁且高度信息化的可视化方法，被实现在多个R图书馆中，比如GGally包中的ggpairs()。相关矩阵图的上对角线显示我们数据中数值变量之间的双变量相关性，或者皮尔逊相关系数。在下对角线中，我们看到了拟合到我们数据的散点图和回归曲线。最后，在主对角线上，我们有每个变量的直方图和密度曲线。
下面的代码是使用Seaborn的PairGrid()绘图和另一个函数进行适应的实现。
如预期的那样，我们的站点之间存在高度相关性，特别是距离较近的站点，比如北温哥华的两个站点。
需要注意的是，为了缩短计算时间，我们的数据被聚合为6小时期间。如果您使用更大的聚合周期来实验这个图表，您将会看到相关系数的增加，因为均值聚合倾向于平滑掉数据中的异常值。
如果您已经了解了时间序列分析，您现在可能会考虑其他值得检查的相关性类型。但首先，我们需要测试我们的时间序列是否具有平稳性。
平稳性
平稳时间序列是指其统计性质随时间不变的时间序列。换句话说，它具有恒定的均值、方差和自相关性，与时间无关[4]。
一些预测模型依赖于时间序列的平稳性，因此在这个探索阶段测试平稳性的重要性就显而易见了。我们的下一个函数将利用statsmodels实现两个常用的平稳性测试：增强迪基-福勒（”ADF”）测试和Kwiatkowski-Phillips-Schmidt-Shin（”KPSS”）测试。
下面是这两个测试的假设。请注意，它们具有相反的零假设，因此我们将创建一个“决策”列，以便更容易地解释它们的结果。您可以在statsmodels文档中详细了解这两个测试。
增强迪基-福勒（ADF）测试假设：
• H0：时间序列样本中存在单位根（非平稳）
• Ha：时间序列样本中不存在单位根（平稳）
Kwiatkowski-Phillips-Schmidt-Shin（KPSS）测试假设：
• H0：数据围绕常数平稳（平稳）
• Ha：时间序列样本中存在单位根（非平稳）
现在我们有一个机会性的问题是应该在哪个尺度上检查平稳性。答案将高度依赖于您如何对数据建模，而全面的探索性分析的目标之一就是帮助您做出这个决策。
为了说明的目的，在以下示例中，我们将查看温哥华国际机场站2016年1月和2022年1月的数据，看看从2016年到2022年之间的行为是否发生了变化。
您可能还记得在我们的缺失值部分中，我们可以使用Pandas的ffill()、bfill()和interpolate()方法来快速插入系列中的中断部分。您可以看到，我在我们的函数中定义了一个专门的fillna参数，以选择这些方法中的任何一个来快速处理缺失值，因为这两个测试只接受完整样本。
现在回到我们的结果。
我们可以看到，对于2016年，两个测试都表明了非稳定性，但对于2022年，结果出现了分歧。《statsmodels文档》清楚地列出了当ADF和KPSS测试一起进行时的结果解释[6]：
• 案例1：两个测试均得出结论，该系列不稳定——该系列不稳定
• 案例2：两个测试均得出结论，该系列是稳定的——该系列是稳定的
• 案例3：KPSS表明是稳定的，ADF表明是非稳定的——该系列是趋势稳定的。需要去除趋势使系列成为严格的稳定系列。对去趋势后的系列进行稳定性检验。
• 案例4：KPSS表明是非稳定的，ADF表明是稳定的——该系列是差分稳定的。需要使用差分使系列成为稳定的。对差分后的系列进行稳定性检验。
如果你对所有4个站点在多个月份上重复这个操作，你会发现案例4在数据中占主导地位。这将引导我们进入关于一阶差分使数据稳定的下一部分。
一阶差分
作为最常见的转换技术之一，对时间序列应用一阶或二阶差分是广泛用于使数据适用于只能用于稳定时间序列的统计模型的方法。在这里，我们将查看在2016年1月期间对先前的示例之一应用的技术。但首先，让我们用我们的plot_sequence()函数查看转换之前的原始数据。
我们可以看到，在这段时间内，期间变化的方差显著变大。平均PM2.5值似乎从较高的值变为较低且更稳定的值。这些是确认该系列非稳定性的一些特征。
同样，Pandas有一个非常方便的方法来对我们的数据进行区分。我们将调用.diff()到我们的DataFrame，并将其实例化为我们数据的一阶差分版本。所以让我们再次绘制同一段时间。
除了仍然波动的方差，现在的数据明显更稳定地围绕着一个平均值。我们可以再次调用我们的stationarity_test()函数来检查差分后的数据的稳定性。
我们看到了这个。我们在全面的探索性时间序列分析中又进行了另一次检查，因为我们现在确认了：
- 我们正在处理非稳定的时间序列。
- 一阶差分是使其稳定的适当的转换技术。
这最后导致我们到了最后一部分。
自相关
一旦我们的数据是平稳的，我们可以研究其他关键时间序列属性：偏自相关和自相关。在正式的术语中：
自相关函数（ACF）度量时间序列滞后值之间的线性关系。换句话说，它度量时间序列与自己的相关性。[2]
偏自相关函数（PACF）度量时间序列中滞后值之间的相关性，但是在移除了中间相关滞后值的影响时。这些被称为混淆因子。[3]
这两个指标可以用称为相关图的统计图可视化。但首先，重要的是对它们有更好的理解。
由于本文侧重于探索性分析，并且这些概念对于统计预测模型非常重要，我将简要介绍解释，但请记住，当处理时间序列时，这些都是非常重要的概念以建立扎实的直觉。
综合阅读，我推荐Kaggle Notebooks大师Leonie Monigatti的优秀文章“时间序列：解释ACF和PACF”。
如上所述，自相关度量时间序列与以前q滞后的自身相关程度。您可以将它看作是数据子集与其自身的副本向后移动q个时期的线性关系的度量。自相关，或者ACF，是确定移动平均（MA）模型的阶数q的重要指标。
另一方面，偏自相关是时间序列与其p滞后版本的相关性，但现在只涉及其直接影响。例如，如果我想检查从t-3到t-1时间段与当前t0值的偏自相关性，我不关心t-3如何影响t-2和t-1，或者t-2如何影响t-1。我将专注于t-3、t-2和t-1对当前时间戳t0的直接影响。偏自相关，或者PACF，是确定自回归（AR）模型的阶数p的重要指标。
现在我们清楚了这些概念，我们可以回到我们的数据。由于这两个指标经常一起分析，我们的最后一个函数将在网格图中结合PACF和ACF图，返回多个变量的相关图。它将使用statsmodels的plot_pacf()和plot_acf()函数，并将它们映射到Matplotlib的subplots()网格。
请注意，除了plot_pacf()图独有的method参数外，statsmodels的这两个函数使用相同的参数。
现在您可以尝试使用数据的不同汇总方式，但请记住，当重新采样时间序列时，每个滞后将代表不同的时间跳跃。出于说明目的，让我们分析2016年1月四个站点的PACF和ACF，使用6小时的汇总数据集。
相关图返回从-1.0到1.0的相关系数以及表示显著性阈值的阴影区域。超出该范围的任何值都应视为具有统计意义。
根据上面的结果，最后我们可以得出以下结论：在6小时汇总中：
- 滞后1、2、3（t-6h、t-12h和t-18h），有时4（t-24h）具有显著的PACF。
- 滞后1和4（t-6h和t-24h）在大多数情况下显示出显著的ACF。
并注意一些最终的良好实践：
- 应避免为具有高精度的大时间序列绘制相关图（例如，为具有小时测量数据的数据集绘制整年的相关图），因为随着样本大小的增加，显著性阈值会趋近于零。
- 我在我们的函数中定义了一个
x_label参数，以便于用每个滞后期所表示的时间段对X轴进行注释。在不包含该信息的相关图中很常见，但是通过轻松访问它可以避免对结果的错误解释。 - Statsmodels 的
plot_acf()和plot_pacf()默认值设置为在图中包含0滞后期的相关系数。由于任何数与自身的相关性总是为1，我将我们的图形设置为从第一个滞后期开始，参数为zero=False。它还改善了Y轴的刻度，使我们实际需要分析的滞后期更容易阅读。
有了这些，我们已经彻底探索了我们的时间序列。通过可视化工具和分析函数的一套工具，我们可以对数据有一个全面的了解。您还学到了探索时间序列数据集的最佳实践以及如何通过高质量的图表简洁而精美地呈现它们。
喜欢这个故事吗？
您可以在VoAGI上关注我，以获取更多关于数据科学、机器学习、可视化和数据分析的文章。
您也可以在LinkedIn上找到我，我在那里分享这些内容的较短版本。
参考资料
[1] “卫生部 – 微粒子（PM 2.5）问题和答案。” 2022年10月14日访问。 https://www.health.ny.gov/environmental/indoors/air/pmq_a.htm.
[2] Peixeiro，Marco。“3.进行随机游走。” 文章。《Python时间序列预测》，30–58。 O’Reilly Media，2022。
[3] Peixeiro，Marco。“5.建模自回归过程。” 文章。《Python时间序列预测》，81–100。 O’Reilly Media，2022。
[4] Peixeiro，Marco。“8.考虑季节性。” 文章。《Python时间序列预测》，156–79。 O’Reilly Media，2022。
[5] 服务，市民部。 “卑诗省数据目录。” 卑诗省。 卑诗省，2022年2月2日。 https://www2.gov.bc.ca/gov/content/data/bc-data-catalogue.
[6] “平稳性和去趋势（ADF/KPSS）。 ”statsmodels。 2022年10月17日访问。 https://www.statsmodels.org/dev/examples/notebooks/generated/stationarity_detrending_adf_kpss.html.