python时间序列特征提取_干货undefined|undefined20个Python教程，掌握时间序列的特征分析（附代码）...

作者 | Selva Prabhakaran

译者 | Tianyu

编辑 | Jane

来源 | AI科技大本营(ID: rgznai100)

【导语】时间序列是指以固定时间为间隔的序列值。本篇教程将教大家用 Python 对时间序列进行特征分析。

1、什么是时间序列？

时间序列是指以固定时间为间隔的、由所观察的值组成的序列。根据观测值的不同频率，可将时间序列分成小时、天、星期、月份、季度和年等时间形式的序列。有时候，你也可以将秒钟和分钟作为时间序列的间隔，如每分钟的点击次数和访客数等等。

为什么我们要对时间序列进行分析呢？

因为当你想对一个序列进行预测时，首先要完成分析这个步骤。除此之外，时间序列的预测也具有极大商业价值，如企业的供求量、网站的访客量以及股票价格等，都是极其重要的时间序列数据。

那么，时间序列分析都包括哪些内容呢？

要做好时间序列分析，必须要理解序列的内在属性，这样才能做出更有意义且精准的预测。

2、如何在 Python 中引入时间序列？

关于时间序列的数据大都存储在 csv 文件或其他形式的表格文件里，且都包含两个列：日期和观测值。

首先我们来看 panda 包里面的 read_csv() 函数，它可以将时间序列数据集(关于澳大利亚药物销售的 csv 文件)读取为 pandas 数据框。增加一个 parse_dates=['date'] 字段，可以把包含日期的数据列解析为日期字段。from dateutil.parser import parseimport matplotlib as mplimport matplotlib.pyplot as pltimport seaborn as snsimport numpy as npimport pandas as pdplt.rcParams.update({'figure.figsize': (10, 7), 'figure.dpi': 120})# Import as Dataframedf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])df.head()

pandas 序列

注意，在 pandas 序列中，'value' 列的位置高于 'date' 列，这表明它是一个 pandas 序列而非数据框。

3、什么是面板数据？

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

不同之处是，除了时间序列，面板数据还包括一个或多个相关变量，这些变量也是在同个时间段内测得的。

面板数据中的列包括有助于预测 y 值的解释变量，这些特征列可用于之后的预测。以下是关于面板数据的例子：# dataset source: https://github.com/rouseguydf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/MarketArrivals.csv')df = df.loc[df.market=='MUMBAI', :]df.head()

时间序列可视化

由于所有的值都为正数，无论使用 y 轴的哪一侧都可以看到增长的趋势。# Import datadf = pd.read_csv('datasets/AirPassengers.csv', parse_dates=['date'])x = df['date'].valuesy1 = df['value'].values# Plotfig, 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()

药品销售的季节性时间序列图

每年的二月份，药品销售会有一次大幅度下跌，三月份会回涨，四月份再次下跌，以此规律循环。很明显，该模式以年为单位，每年循环往复。

不过，随着年份的变化，药品销售总体呈上升趋势。你可以选择使用箱型图将这一趋势进行可视化，可以方便看出每一年的变化。同样地，你也可以按月份绘制箱型图，来观察每个月的变化。

按月份(季节)和年份绘制箱型图：你可以将数据处理成以季节为时间间隔，然后观察特定年份内值的分布，也可以将全部时间的数据进行对比。# Import Datadf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')df.reset_index(inplace=True)# Prepare datadf['year'] = [d.year for d in df.date]df['month'] = [d.strftime('%b') for d in df.date]years = df['year'].unique()# Draw Plotfig, 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 Titleaxes[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

时间序列的模式

另一个需要考虑的方面是周期性模式。当序列中的上升和下降，不是按日历中的特定时间间隔发生时，就会出现这种情况。注意不要把“周期”作用和“季节”作用混淆。

那么，如何区分“周期”和“季节”呢？

如果序列中的模式不是以日历中特定间隔循环出现的，那么就是周期。因为与季节性不同，周期作用通常受到商业或社会经济等因素的影响。

6、加法与乘法时间序列

根据趋势和季节的固有属性，一个时间序列可以被建模为加法模型或乘法模型，也就是说，序列中的值可以用各个成分的加和或乘积来表示：

加法时间序列：

值 = 基准 + 趋势 + 季节 + 残差

乘法时间序列：

值 = 基准 x 趋势 x 季节 x 残差

7、如何将时间序列的成分分解出来？

通过将一个时间序列视为基准、趋势、季节指数及残差的加法或乘法组合，你可以对时间序列进行经典分解。statsmodels 的 seasonal_decompose 函数可以使这一过程非常容易。from statsmodels.tsa.seasonal import seasonal_decomposefrom dateutil.parser import parse# Import Datadf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date')# Multiplicative Decompositionresult_mul = seasonal_decompose(df['value'], model='multiplicative', extrapolate_trend='freq')# Additive Decompositionresult_add = seasonal_decompose(df['value'], model='additive', extrapolate_trend='freq')# Plotplt.rcParams.update({'figure.figsize': (10,10)})result_mul.plot().suptitle('Multiplicative Decompose', fontsize=22)result_add.plot().suptitle('Additive Decompose', fontsize=22)plt.sho

平稳和非平稳时间序列

上图来自 R 语言的 TSTutorial 教程。

为什么说序列平稳很重要呢？

我会对此稍微做一些解释，但要清楚一点，通过采用合适的变换，我们近乎可以将任何时间序列变得平稳。大多数统计预测方法最初都是为平稳时间序列而设计的。预测过程的第一步是通过一些变换，来将非平稳序列变成平稳序列。

9、如何将时间序列变平稳？

你可以通过以下几种方式得到平稳序列：求序列的差分

求序列的 log 值

求序列的 n 次方根

把上面三种方法相结合

将时间序列平稳化最普遍且便捷的方法是对序列进行差分运算，至少执行一次，直到序列趋于平稳。

那么什么是差分呢？

如果 Y_t 为时间 't' 对应的值，那么第一个差分值为 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]

10、为什么要在预测前将非平稳序列变平稳？

对平稳序列进行预测要相对容易一些，预测结果也更可信。其中一个重要原因是，自回归预测模型本质上是线性回归模型，将序列自身的滞后作为预测因子。

如果预测因子之间互不相关，线性回归的效果最优。那么将序列平稳化就可以解决这一问题，因为它可以去除任何持久的自相关性，所以可以使预测模型中的预测因子近乎独立。

现在我们知道了使序列平稳化的重要性，那么应该如何检查一个序列是否平稳呢？

11、如何测试平稳性？

我们可以像之前那样，通过绘制序列图来观察一个序列是否平稳。

另一种方法是将序列分解成两个或多个连续的部分，并求其统计值，如平均值、方差和自相关系数。如果这些统计值间的差异很大，那么该序列大概率不是平稳序列。

尽管如此，你仍需要一种方法来定量地判断某个序列是否平稳。一个名为“单位根检验”的统计检验方法可以解决这一问题。

单位根检验有如下几种方法：ADF Test (Augmented Dickey Fuller test)

KPSS Test (Kwiatkowski-Phillips-Schmidt-Shin) — 趋势平稳

PP Test (Philips Perron test)

最常用的方法是 ADF test，零假设为：用时间序列对单位根进行处理，结果是不平稳。如果 P 值小于显著性水平 0.05，则拒绝零假设，即不成立。

另一方面，KPSS test 可用来检测趋势平稳性。零假设与 P 值含义都与 ADH test 相反。以下代码基于 Python 的 statsmodels 包执行了两种检测方法：from statsmodels.tsa.stattools import adfuller, kpssdf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])# ADF Testresult = 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 Testresult = 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}'from statsmodels.tsa.stattools import adfuller, kpssdf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])# ADF Testresult = 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 Testresult = 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}'

12、白噪声与平稳序列的差别是什么？

和平稳序列一样，白噪声也是关于时间的函数，它的均值和方差始终不变。但区别在于，白噪声是完全随机的，且均值为零。

白噪声没有任何模式。如果你把调频广播的声波讯号想象成时间序列，调频道时的空白声音就是白噪声。

从数学上来讲，一个完全随机且均值为零的序列就是白噪声。randvals = np.random.randn(1000)pd.Series(randvals).plot(title='Random White Noise', color='k')

通过减掉最小二乘拟合线对时间序列去趋势# Using statmodels: Subtracting the Trend Component.from statsmodels.tsa.seasonal import seasonal_decomposedf = 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.trendplt.plot(detrended)plt.title('Drug Sales detrended by subtracting the trend component', fontsize=16)

对时间序列去季节

15、如何检测时间序列的季节性？

一般方法是画出序列图，观察固定的时间间隔是否有重复模式出现。因此，季节性的类型由时钟或日历决定：一天中的小时

月份中的日期

星期

月份

年份

不过，如果你想对季节性做一个明确的检验，可以使用自相关函数 (ACF) 图，接下来的部分会做相关详细介绍。当季节模式明显时，ACF 图中季节窗口的整数倍处会反复出现特定的尖峰。

例如，药品销售的时间序列是月份序列，每年会出现重复的模式，你会在第 12、24、36 个序列值处看到尖峰。

必须要提醒你的是，现实生活中的数据集很少会出现特别明显的模式，可能会被一些噪声污染，所以需要更加仔细观察其中的模式。from pandas.plotting import autocorrelation_plotdf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv')# Draw Plotplt.rcParams.update({'figure.figsize':(9,5), 'figure.dpi':120})autocorrelation_plot(df.value.tolist())

# # Generate datasetfrom scipy.interpolate import interp1dfrom sklearn.metrics import mean_squared_errordf_orig = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'], index_col='date').head(100)df = pd.read_csv('datasets/a10_missings.csv', parse_dates=['date'], index_col='date')fig, axes = plt.subplots(7, 1, sharex=True, figsize=(10, 12))plt.rcParams.update({'xtick.bottom' : False})## 1. Actual -------------------------------df_orig.plot(title='Actual', ax=axes[0], label='Actual', color='red', style=".-")df.plot(title='Actual', ax=axes[0], label='Actual', color='green', style=".-")axes[0].legend(["Missing Data", "Available Data"])## 2. Forward Fill --------------------------df_ffill = df.ffill()error = np.round(mean_squared_error(df_orig['value'], df_ffill['value']), 2)df_ffill['value'].plot(title='Forward Fill (MSE: ' + str(error) +")", ax=axes[1], label='Forward Fill', style=".-")## 3. Backward Fill -------------------------df_bfill = df.bfill()error = np.round(mean_squared_error(df_orig['value'], df_bfill['value']), 2)df_bfill['value'].plot(title="Backward Fill (MSE: " + str(error) +")", ax=axes[2], label='Back Fill', color='firebrick', style=".-")## 4. Linear Interpolation ------------------df['rownum'] = np.arange(df.shape[0])df_nona = df.dropna(subset = ['value'])f = interp1d(df_nona['rownum'], df_nona['value'])df['linear_fill'] = f(df['rownum'])error = np.round(mean_squared_error(df_orig['value'], df['linear_fill']), 2)df['linear_fill'].plot(title="Linear Fill (MSE: " + str(error) +")", ax=axes[3], label='Cubic Fill', color='brown', style=".-")## 5. Cubic Interpolation --------------------f2 = interp1d(df_nona['rownum'], df_nona['value'], kind='cubic')df['cubic_fill'] = f2(df['rownum'])error = np.round(mean_squared_error(df_orig['value'], df['cubic_fill']), 2)df['cubic_fill'].plot(title="Cubic Fill (MSE: " + str(error) +")", ax=axes[4], label='Cubic Fill', color='red', style=".-")# Interpolation References:# https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html# https://docs.scipy.org/doc/scipy/reference/interpolate.html## 6. Mean of 'n' Nearest Past Neighbors ------def knn_mean(ts, n):out = np.copy(ts)for i, val in enumerate(ts):if np.isnan(val):n_by_2 = np.ceil(n/2)lower = np.max([0, int(i-n_by_2)])upper = np.min([len(ts)+1, int(i+n_by_2)])ts_near = np.concatenate([ts[lower:i], ts[i:upper]])out[i] = np.nanmean(ts_near)return outdf['knn_mean'] = knn_mean(df.value.values, 8)error = np.round(mean_squared_error(df_orig['value'], df['knn_mean']), 2)df['knn_mean'].plot(title="KNN Mean (MSE: " + str(error) +")", ax=axes[5], label='KNN Mean', color='tomato', alpha=0.5, style=".-")## 7. Seasonal Mean ----------------------------def seasonal_mean(ts, n, lr=0.7):"""Compute the mean of corresponding seasonal periodsts: 1D array-like of the time seriesn: Seasonal window length of the time series"""out = np.copy(ts)for i, val in enumerate(ts):if np.isnan(val):ts_seas = ts[i-1::-n] # previous seasons onlyif np.isnan(np.nanmean(ts_seas)):ts_seas = np.concatenate([ts[i-1::-n], ts[i::n]]) # previous and forwardout[i] = np.nanmean(ts_seas) * lrreturn outdf['seasonal_mean'] = seasonal_mean(df.value, n=12, lr=1.25)error = np.round(mean_squared_error(df_orig['value'], df['seasonal_mean']), 2)df['seasonal_mean'].plot(title="Seasonal Mean (MSE: " + str(error) +")", ax=axes[6], label='Seasonal Mean', color='blue', alpha=0.5, s

ACF 和 PACF

18、如何计算偏自相关系数？

序列滞后 k 处的偏自相关系数是 Y 的自回归方程的滞后系数。Y 的自回归方程是指 Y 以自己的滞后作为预测因子的线性回归。

举个例子，如果 Y_t 为当前序列，Y_t-1 即为滞后期为 1 的 Y 值，那么滞后期为 3 处 (Y_t-3) 的偏自相关系数是下面方程中的 α3：

药品销售滞后图

时间序列的平滑处理

22、如何用格兰杰因果关系检验判断一个时间序列是否可以预测另一个？

格兰杰因果关系检验可用作检测一个时间序列是否可以用来预测另一个序列。那么格兰杰因果关系检验是如何进行运算的呢？

该检验基于一个假设，即 X 导致了 Y 的发生，那么基于 Y 的先前值与 X 的先前值得到的 Y 的预测值，要优于仅根据 Y 本身得到的预测值。

statsmodels 包提供了便捷的检验方法。它可以接受一个二维数组，其中第一列为值，第二列为预测因子。

零假设为：第二列的序列与第一列不存在格兰杰因果关系。如果 P 值小于显著性水平 0.05，就可以拒绝零假设，从而知道 X 的滞后是起作用的。

第二个参数 maxlag 表示该检验中涉及了 Y 的多少个滞后期。from statsmodels.tsa.stattools import grangercausalitytestsdf = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])df['month'] = df.date.dt.monthgrangercausalitytests(df[['value', 'month']], maxlag=2)Granger Causalitynumber of lags (no zero) 1ssr based F test: F=54.7797 , p=0.0000 , df_denom=200, df_num=1ssr based chi2 test: chi2=55.6014 , p=0.0000 , df=1likelihood ratio test: chi2=49.1426 , p=0.0000 , df=1parameter F test: F=54.7797 , p=0.0000 , df_denom=200, df_num=1Granger Causalitynumber of lags (no zero) 2ssr based F test: F=162.6989, p=0.0000 , df_denom=197, df_num=2ssr based chi2 test: chi2=333.6567, p=0.0000 , df=2likelihood ratio test: chi2=196.9956, p=0.0000 , df=2parameter F test:         F=162.6989, p=0.0000  , df_denom=197, df_num=2

在上面的例子中，全部测试结果中的 P 值都为零，说明 'month' 可用作预测航班的乘客数量。原文链接：

https://www.machinelearningplus.com/time-series/time-series-analysis-python/

(*本文为 Python大本营转载文章，转载请联系1092722531)