推荐收藏时间序列分析全面指南(附Python代码)

Posted 我爱Python数据挖掘

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了推荐收藏时间序列分析全面指南(附Python代码)相关的知识,希望对你有一定的参考价值。

大家好,时间序列是在规律性时间间隔上记录的观测值序列。本文我将带你了解在 Python 中分析给定时间序列的特征的21个全过程。内容较长,建议收藏、点赞、关注。

内容

本文部分技术来自技术群小伙伴的分享,加入按照如下方式

目前开通了技术交流群,群友已超过3000人,添加时最好的备注方式为:来源+兴趣方向,方便找到志同道合的朋友
方式①、添加微信号:dkl88191,备注:来自CSDN+技术交流
方式②、微信搜索公众号:Python学习与数据挖掘,后台回复:加群+CSDN

1. 什么是时间序列?

时间序列是在规律性时间间隔记录的观测值序列。

依赖于观测值的频率,典型的时间序列可分为每小时、每天、每周、每月、每季度和每年为单位记录。有时,你可能也会用到以秒或者分钟为单位的时间序列,比如,每分钟用户点击量和访问量等等。

1.1 为什么要分析时间序列呢?

因为它是你做序列预测前的一步准备过程。而且,时间序列预测拥有巨大的商业重要性,因为对商业来说非常重要的需求和销量、网站访问人数、股价等都是时间序列数据。

1.2 所以时间序列分析包括什么内容呢?

时间序列分析包括理解序列内在本质的多个方面以便于你可更好地了解如何做出有意义并且精确的预测。

2. 如何在Python中导入时间序列?


所以怎样导入时间序列数据呢?典型的时间序列数据以.csv格式或者其他表格形式存储,包括两列:日期和测量值。

让我们用pandas包里的read.csv()读取时间序列数据(一个澳大利亚药品销售的csv文件)作为一个pandas数据框。加入parse_dates=[‘date’]参数将会把日期列解析为日期字段。

from dateutil.parser import parse
import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import pandas as pd
plt.rcParams.update('figure.figsize': (10, 7), 'figure.dpi': 120)

# Import as Dataframe
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])
df.head()

数据框时间序列

此外,你也可以将其导入为date作为索引的pandas序列。你只需要固定pd.read_csv()里的index_col参数。

ser = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')
ser.head()

时间序列

注意,在此序列当中,‘value’列的位置高于date以表明它是一个序列。

3. 什么是面板数据?

面板数据也是基于时间的数据集。

差异在于,除了时间序列,它也包括同时测量的一个或多个相关变量。

通常来看,面板数据当中的列包括了有助于预测Y的解释型变量,假设这些列将在未来预测阶段有用。

面板数据的例子如下:

# dataset source: https://github.com/rouseguy
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/MarketArrivals.csv')
df = df.loc[df.market=='MUMBAI', :]
df.head()

面板数据

4. 时间序列可视化


让我们用matplotlib来对序列进行可视化。

# Time series data source: fpp pacakge in R.
import matplotlib.pyplot as plt
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')

# Draw Plot
def plot_df(df, x, y, title="", xlabel='Date', ylabel='Value', dpi=100):
    plt.figure(figsize=(16,5), dpi=dpi)
    plt.plot(x, y, color='tab:red')
    plt.gca().set(title=title, xlabel=xlabel, ylabel=ylabel)
    plt.show()

plot_df(df, x=df.index, y=df.value, title='Monthly anti-diabetic drug sales in  Australia from 1992 to 2008.')

时间序列可视化

因为所有的值都是正值,你可以在Y轴的两侧进行显示此值以强调增长。

# Import data
df = pd.read_csv('datasets/AirPassengers.csv', parse_dates=['date'])
x = df['date'].values
y1 = df['value'].values

# Plot
fig, ax = plt.subplots(1, 1, figsize=(16,5), dpi= 120)
plt.fill_between(x, y1=y1, y2=-y1, alpha=0.5, linewidth=2, color='seagreen')
plt.ylim(-800, 800)
plt.title('Air Passengers (Two Side View)', fontsize=16)
plt.hlines(y=0, xmin=np.min(df.date), xmax=np.max(df.date), linewidth=.5)
plt.show()

航空乘客数据——两侧序列

因为这是一个月度时间序列,每年遵循特定的重复模式,你可以把每年作为一个单独的线画在同一张图上。这可以让你同时比较不同年份的模式。

4.1 时间序列的季节图

# Import Data
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')
df.reset_index(inplace=True)

# Prepare data
df['year'] = [d.year for d in df.date]
df['month'] = [d.strftime('%b') for d in df.date]
years = df['year'].unique()

# Prep Colors
np.random.seed(100)
mycolors = np.random.choice(list(mpl.colors.XKCD_COLORS.keys()), len(years), replace=False)

# Draw Plot
plt.figure(figsize=(16,12), dpi= 80)
for i, y in enumerate(years):
    if i > 0:        
        plt.plot('month', 'value', data=df.loc[df.year==y, :], color=mycolors[i], label=y)
        plt.text(df.loc[df.year==y, :].shape[0]-.9, df.loc[df.year==y, 'value'][-1:].values[0], y, fontsize=12, color=mycolors[i])

# Decoration
plt.gca().set(xlim=(-0.3, 11), ylim=(2, 30), ylabel='$Drug Sales$', xlabel='$Month$')
plt.yticks(fontsize=12, alpha=.7)
plt.title("Seasonal Plot of Drug Sales Time Series", fontsize=20)
plt.show()

药品销售的季节图

每年二月会迎来药品销售的急速下降,而在三月会再度上升,接下来的4月又开始下降,以此类推。很明显,该模式在特定的某一年中重复,且年年如此。

然而,随着年份推移,药品销售整体增加。你可以很好地看到该趋势并且在年份箱线图当中看到它是怎样变化的。同样地,你也可以做一个月份箱线图来可视化月度分布情况。

4.2 月度(季节性)箱线图和年度(趋势)分布

你可以季节间隔将数据分组,并看看在给定的年份或月份当中值是如何分布的,以及随时间推移它们是如何比较的。

# Import Data
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')
df.reset_index(inplace=True)

# Prepare data
df['year'] = [d.year for d in df.date]
df['month'] = [d.strftime('%b') for d in df.date]
years = df['year'].unique()

# Draw Plot
fig, axes = plt.subplots(1, 2, figsize=(20,7), dpi= 80)
sns.boxplot(x='year', y='value', data=df, ax=axes[0])
sns.boxplot(x='month', y='value', data=df.loc[~df.year.isin([1991, 2008]), :])

# Set Title
axes[0].set_title('Year-wise Box Plot\\n(The Trend)', fontsize=18); 
axes[1].set_title('Month-wise Box Plot\\n(The Seasonality)', fontsize=18)
plt.show()

年度和月度箱线图

箱线图将年度和月度分布变得很清晰。并且,在阅读箱线图当中,12月和1月明显有更高的药品销售量,可被归因于假期折扣季。

到目前为止,我们已经看到了识别模式的相似之处。现在怎样才能从通常模式当中找到离群值呢?

5. 时间序列的模式


任何时间序列都可以被分解为如下的部分:基线水平+趋势+季节性+误差

当在时间序列当中观测到增加或降低的斜率时,即可观测到相应的趋势。然而季节性只有在由于季节性因素导致不同的重复模式在规律性的间隔之间被观测到时才能发现。可能是由于当年的特定月份,特定月份的某一天、工作日或者甚至是当天某个时间。

然而,并不是所有时间序列必须有一个趋势和/或季节性。时间序列可能没有不同的趋势但是有一个季节性。反之亦然。

所以时间序列可以被看做是趋势、季节性和误差项的整合。

fig, axes = plt.subplots(1,3, figsize=(20,4), dpi=100)
pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/guinearice.csv', parse_dates=['date'], index_col='date').plot(title='Trend Only', legend=False, ax=axes[0])

pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/sunspotarea.csv', parse_dates=['date'], index_col='date').plot(title='Seasonality Only', legend=False, ax=axes[1])

pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/AirPassengers.csv', parse_dates=['date'], index_col='date').plot(title='Trend and Seasonality', legend=False, ax=axes[2])

时间序列中的模式

另一个需要考虑的方面是循环的行为。当序列当中上升和下降模式并不在固定的日历间隔出现时,就会出现循环的行为。需注意不要混淆循环的效应和季节的效应。

所以,怎样区分循环的和季节性的模式呢?

如果模式不是基于固定的日历频率,那它就是循环的。因为,循环效应不像季节性那样受到商业和其他社会经济因素的影响。

6. 时间序列的加法和乘法


基于趋势和季节性的本质,时间序列以加法或乘法的形式建模,其中序列里的每个观测值可被表达为成分的和或者积:

**加法时间序列:**值=基线水平+趋势+季节性+误差

**乘法时间序列:**值=基线水平*趋势*季节性*误差

7. 怎样分解时间序列的成分?


你可以通过将序列作基线水平,趋势,季节性指数和残差的加法或乘法组合来实现一个经典的时间序列分解。

statsmodels包里的seasonal_decompose使用起来非常方便。

from statsmodels.tsa.seasonal import seasonal_decompose
from dateutil.parser import parse

# Import Data
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')

# Multiplicative Decomposition
result_mul = seasonal_decompose(df['value'], model='multiplicative', extrapolate_trend='freq')

# Additive Decomposition
result_add = seasonal_decompose(df['value'], model='additive', extrapolate_trend='freq')

# Plot
plt.rcParams.update('figure.figsize': (10,10))
result_mul.plot().suptitle('Multiplicative Decompose', fontsize=22)
result_add.plot().suptitle('Additive Decompose', fontsize=22)
plt.show()

加法和乘法分解

在序列开始时,设置extrapolate_trend=‘freq’ 来注意趋势和残差中缺失的任何值。

如果你仔细看加法分解当中的残差,它有一些遗留模式。乘法分解看起来非常随意,这很好。所以理想状况下,乘法分解应该在这种特定的序列当中优先选择。

趋势,季节性和残差成分的数值输出被存储在result_mul 当中。让我们提取它们并导入数据框中。

# Extract the Components ----
# Actual Values = Product of (Seasonal * Trend * Resid)
df_reconstructed = pd.concat([result_mul.seasonal, result_mul.trend, result_mul.resid, result_mul.observed], axis=1)
df_reconstructed.columns = ['seas', 'trend', 'resid', 'actual_values']
df_reconstructed.head()

如果你检查一下seas, trend 和 resid列的乘积,应该确实等于actual_values。

8. 平稳和非平稳时间序列


平稳性是时间序列的属性之一。平稳序列的值不是时间的函数。

也就是说,这种序列的统计属性例如均值,方差和自相关是随时间不变的常数。序列的自相关只是与前置值的相关,之后会详细介绍。

平稳时间序列也没有季节效应。

所以如何识别一个序列是否平稳呢?让我们通过实例来展示一下:

平稳和非平稳时间序列

上图来自R语言的 TSTutorial。

所以为什么平稳序列是重要的呢?为什么我要提到它?

我将展开讲一下,但是要理解它只是有可能通过使用特定的转换方法实现任何时间序列的平稳化。大多数统计预测方法都用于平稳时间序列。预测的第一步通常是做一些转换将非平稳数据转化为平稳数据。

9. 如何获取平稳的时间序列?


你可以通过以下步骤实现序列的平稳化:

1. 差分序列(一次或多次);

2. 对序列值进行log转换;

3. 对序列值去n次根式值;

4. 结合上述方法。

实现数据平稳化最常见也最方便的方法是对序列进行差分至少一次,直到它变得差不多平稳为止。

9.1 所以什么是差分?

如果Y_t是t时刻的Y值,那么第一次差分Y = Yt – Yt-1。在简化的格式当中,差分序列就是从当前值中减去下一个值。

如果第一次差分不能使数据平稳,你可以第二次差分,以此类推。

例如,考虑如下序列: [1, 5, 2, 12, 20]

一次差分: [5-1, 2-5, 12-2, 20-12] = [4, -3, 10, 8]

二次差分: [-3-4, -10-3, 8-10] = [-7, -13, -2]


9.2为什么要在预测之前将非平稳数据平稳化?

预测平稳序列相对容易,预测也相对更可靠。

一个重要的原因是自回归预测模型必须是利用序列自身的滞后量作为预测变量的线性回归模型。

我们知道线性回归在预测变量(X变量)与其他变量不相关时效果最佳。所以序列平稳化也因为移除所有持续的自相关而解决了这个问题,因此使得模型中的预测变量(序列的滞后值)几乎独立。

现在我们已经建立了序列平稳化非常重要的概念,那怎样检验给定序列是否平稳化呢?

10. 怎样检验平稳性?


序列的平稳性可以通过之前我们提到的序列图看出来。

另外一种方法是将序列分成2或多个连续的部分,计算概要统计量例如均值,方差和自相关。如果统计量显著差异,序列可能不是平稳的。

尽管如此,你需要一个方法来从量化的角度判断一个给定序列是否平稳。可以通过‘Unit Root Tests单位根检验’来实现。这里有多种变式,但这些检验都是用来检测时间序列是否非平稳并且拥有一个单位根。

有多种单位根检验的具体应用:

  1. 增广迪基·富勒检验(ADF Test);

2. 科维亚特夫斯基-菲利普斯-施密特-辛-KPSS检验(趋势平稳性);

3. 菲利普斯 佩龙检验(PP Test)。

最常用的是ADF检验,零假设是时间序列只有一个单位根并且非平稳。所以ADF检验p值小于0.05的显著性水平,你拒绝零假设。

KPSS检验,另一方面,用于检验趋势平稳性。零假设和p值解释与ADH检验相反。下面的代码使用了python中的statsmodels包来做这两种检验。

from statsmodels.tsa.stattools import adfuller, kpss
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])

# ADF Test
result = adfuller(df.value.values, autolag='AIC')
print(f'ADF Statistic: result[0]')
print(f'p-value: result[1]')
for key, value in result[4].items():
    print('Critial Values:')
    print(f'   key, value')

# KPSS Test
result = kpss(df.value.values, regression='c')
print('\\nKPSS Statistic: %f' % result[0])
print('p-value: %f' % result[1])
for key, value in result[3].items():
    print('Critial Values:')
    print(f'   key, value')
ADF Statistic: 3.14518568930674
p-value: 1.0
Critial Values:
   1%, -3.465620397124192
Critial Values:
   5%, -2.8770397560752436
Critial Values:
   10%, -2.5750324547306476

KPSS Statistic: 1.313675
p-value: 0.010000
Critial Values:
   10%, 0.347
Critial Values:
   5%, 0.463
Critial Values:
   2.5%, 0.574
Critial Values:
   1%, 0.739

11. 白噪音和平稳序列的差异是什么?


如平稳序列,白噪音也不是时间的函数,它的均值和方差并不随时间变化。但是它与平稳序列的差异在于,白噪音完全随机,均值为0。

无论怎样,在白噪音当中是没有特定模式的。如果你将FM广播的声音信号作为时间序列,你在频道之间的频段听到的空白声就是白噪音。

从数学上来看,均值为0的完全随机的数字序列是白噪音。

randvals = np.random.randn(1000)
pd.Series(randvals).plot(title='Random White Noise', color='k')

随机白噪音

12. 怎样将时间序列去趋势化?


对时间序列去趋势就是从时间序列当中移除趋势成分。但是如何提取趋势呢?有以下几个方法。

1. 从时间序列当中减去最优拟合线。最佳拟合线可从以时间步长为预测变量获得的线性回归模型当中获得。对更复杂的模型,你可以使用模型中的二次项(x^2);

2. 从我们之前提过的时间序列分解当中减掉趋势成分;

3. 减去均值;

4. 应用像Baxter-King过滤器(statsmodels.tsa.filters.bkfilter)或者Hodrick-Prescott 过滤器 (statsmodels.tsa.filters.hpfilter)来去除移动的平均趋势线或者循环成分。

让我们来用一下前两种方法。

# Using scipy: Subtract the line of best fit
from scipy import signal
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])
detrended = signal.detrend(df.value.values)
plt.plot(detrended)
plt.title('Drug Sales detrended by subtracting the least squares fit', fontsize=16)

通过减去最小二乘拟合来对时间序列去趋势化

# Using statmodels: Subtracting the Trend Component.
from statsmodels.tsa.seasonal import seasonal_decompose
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')
result_mul = seasonal_decompose(df['value'], model='multiplicative', extrapolate_trend='freq')
detrended = df.value.values - result_mul.trend
plt.plot(detrended)
plt.title('Drug Sales detrended by subtracting the trend component', fontsize=16)

通过减去趋势成分来去趋势化

13. 怎样对时间序列去季节化?


这里有多种方法对时间序列去季节化。以下就有几个:

1.

以上是关于推荐收藏时间序列分析全面指南(附Python代码)的主要内容,如果未能解决你的问题,请参考以下文章

50个数据可视化最有价值的图表(附完整Python代码,建议收藏)

硬核,288页Python核心知识笔记(附思维导图,强烈推荐收藏)

推荐一本适合初学者全面自学python的书(附电子书)

Hive从入门到精通,HQL硬核整理四万字,全面总结,附详细解析,赶紧收藏吧!!

1.8W字TypeScript入门指南:附大量代码实例(收藏!)

推荐这9个 Python 初学者网站值得收藏