动手实战 | 新拿到一批时序数据可以做哪些分析？

VachelHu

发布于 2023-09-02 11:58:12

2920

发布于 2023-09-02 11:58:12

文章被收录于专栏：时序人

对时间序列数据进行分析在很多工业场景里都能遇到。依赖于观测值的频率，典型的时间序列可分为每小时、每天、每周、每月、每季度和每年为单位记录。有时，你可能也会用到以秒或者分钟为单位的时间序列，比如，每分钟用户点击量和访问量等等。

分析时间序列数据非常重要，因为它是你做序列预测前的必不可少的准备过程。理解序列本质的方方面面可以帮助你更好地了解如何做出有意义并且精确的预测。

那么，具体该如何操作呢？本文将为大家简要介绍。点击原文获取英文博客原址。

如何在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以表明它是一个序列。

什么是面板数据？

面板数据也是基于时间的数据集。差异在于，除了时间序列，它也包括同时测量的一个或多个相关变量。通常来看，面板数据当中的列包括了有助于预测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()

时间序列可视化

让我们用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月明显有更高的药品销售量，可被归因于假期折扣季。

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

时间序列的模式

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

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

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

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

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])

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

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

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

时间序列的加法和乘法

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

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

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

怎样分解时间序列的成分？

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

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。

怎样将时间序列去趋势化？

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

从时间序列当中减去最优拟合线。最佳拟合线可从以时间步长为预测变量获得的线性回归模型当中获得。对更复杂的模型，你可以使用模型中的二次项（x^2）；
从我们之前提过的时间序列分解当中减掉趋势成分；
减去均值；
应用像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)

怎样将时间序列去季节化？

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

取一个以长度为季节窗口的移动平均线。这将在这个过程中使序列变得平滑；
序列季节性差分（从当前值当中减去前一季节的值）；
将序列值除以从STL分解当中获得的季节性指数。

如果除以季节性指数后仍没办法得到良好的结果，再试一下序列对数转换然后再做。你之后可以通过去指数恢复到原始尺度。

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

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

# Deseasonalize
deseasonalized = df.value.values / result_mul.seasonal

# Plot
plt.plot(deseasonalized)
plt.title('Drug Sales Deseasonalized', fontsize=16)
plt.plot()

怎样检验时间序列的季节性？

常见的方法是绘制序列并在固定的时间间隔内检查可重复的模式。所以，季节性的类型由钟表或日历决定：

一天的每个小时
一月的每天
每周
每月
每年
...

然而，如果你想要一个更权威的季节性检验，使用自回归函数（ACF）图。更多关于自回归的信息将在下一部分介绍。但是当强季节性模式出现时，ACF图通常揭示了在季节窗的倍数处明显的重复峰值。

例如，药品销售时间序列是每年都有重复模式的一个月度序列。所以，你可以看到第12,24和36条线等的峰值。

from pandas.plotting import autocorrelation_plot
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv')

# Draw Plot
plt.rcParams.update({'figure.figsize':(9,5), 'figure.dpi':120})
autocorrelation_plot(df.value.tolist())

除此之外，如果你想做统计检验，CHT检验可以检验季节性差异是否对序列平稳化有必要。

怎样估计时间序列的预测能力？

时间序列越有规律性和重复性的模式，越容易被预测。“近似熵”可用于量化时间序列波动的规律性和不可预测性。近似熵越高，预测越难。

另一个更好的选项是“样本熵”。样本熵类似与近似熵，但是在估计小时间序列的复杂性上结果更一致。例如，较少样本点的随机时间序列 “近似熵”可能比一个更规律的时间序列更低，然而更长的时间序列可能会有一个更高的“近似熵”。

样本熵可以很好地处理这个问题。请看如下演示：

# https://en.wikipedia.org/wiki/Approximate_entropy
ss = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/sunspotarea.csv')
a10 = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv')
rand_small = np.random.randint(0, 100, size=36)
rand_big = np.random.randint(0, 100, size=136)

def ApEn(U, m, r):
    """Compute Aproximate entropy"""
    def _maxdist(x_i, x_j):
        return max([abs(ua - va) for ua, va in zip(x_i, x_j)])

    def _phi(m):
        x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)]
        C = [len([1 for x_j in x if _maxdist(x_i, x_j) <= r]) / (N - m + 1.0) for x_i in x]
        return (N - m + 1.0)**(-1) * sum(np.log(C))

    N = len(U)
    return abs(_phi(m+1) - _phi(m))

print(ApEn(ss.value, m=2, r=0.2*np.std(ss.value)))     # 0.651
print(ApEn(a10.value, m=2, r=0.2*np.std(a10.value)))   # 0.537
print(ApEn(rand_small, m=2, r=0.2*np.std(rand_small))) # 0.143
print(ApEn(rand_big, m=2, r=0.2*np.std(rand_big)))     # 0.716
0.6514704970333534
0.5374775224973489
0.0898376940798844
0.7369242960384561
# https://en.wikipedia.org/wiki/Sample_entropy
def SampEn(U, m, r):
    """Compute Sample entropy"""
    def _maxdist(x_i, x_j):
        return max([abs(ua - va) for ua, va in zip(x_i, x_j)])

    def _phi(m):
        x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)]
        C = [len([1 for j in range(len(x)) if i != j and _maxdist(x[i], x[j]) <= r]) for i in range(len(x))]
        return sum(C)

    N = len(U)
    return -np.log(_phi(m+1) / _phi(m))

print(SampEn(ss.value, m=2, r=0.2*np.std(ss.value)))      # 0.78
print(SampEn(a10.value, m=2, r=0.2*np.std(a10.value)))    # 0.41
print(SampEn(rand_small, m=2, r=0.2*np.std(rand_small)))  # 1.79
print(SampEn(rand_big, m=2, r=0.2*np.std(rand_big)))      # 2.42
0.7853311366380039
0.41887013457621214
inf
2.181224235989778

del sys.path[0]

如何得知一个时序有助于预测另一个时序？

Granger因果检验被用于检验是否一个时间序列可以预测另一个序列。Granger因果检验是如何工作的？

它基于如果X引起Y的变化，Y基于之前的Y值和之前的X值的预测效果要优于仅基于之前的Y值的预测效果。所以需要了解Granger因果检验不能应用于Y的滞后量引起Y自身的变化的情况，而通常仅用于外源变量（不是Y的滞后量）。

它在statsmodel包中得到了很好的实现。它采纳2列数据的二维数组作为主要参数，被预测值是第一列，而预测变量（X）在第二列。

零假设检验：第二列的序列不能Granger预测第一列数据。如果p值小于显著性水平（0.05），你可以拒绝零假设并得出结论：X的滞后量确实有用。

第二个参数maxlag决定有多少Y的滞后量应该纳入检验当中。

from statsmodels.tsa.stattools import grangercausalitytests
df = pd.read_csv('https://raw.githubusercontent.com/selva86/datasets/master/a10.csv', parse_dates=['date'])
df['month'] = df.date.dt.month
grangercausalitytests(df[['value', 'month']], maxlag=2)
Granger Causality
number of lags (no zero) 1
ssr based F test:         F=54.7797 , p=0.0000  , df_denom=200, df_num=1
ssr based chi2 test:   chi2=55.6014 , p=0.0000  , df=1
likelihood ratio test: chi2=49.1426 , p=0.0000  , df=1
parameter F test:         F=54.7797 , p=0.0000  , df_denom=200, df_num=1

Granger Causality
number of lags (no zero) 2
ssr based F test:         F=162.6989, p=0.0000  , df_denom=197, df_num=2
ssr based chi2 test:   chi2=333.6567, p=0.0000  , df=2
likelihood ratio test: chi2=196.9956, p=0.0000  , df=2
parameter F test:         F=162.6989, p=0.0000  , df_denom=197, df_num=2

在上述例子当中，实际上所有检验的p值只能无限接近于0。所以“月份”实际上可以用于预测航空乘客的数量。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2023-05-23，如有侵权请联系 cloudcommunity@tencent.com 删除

csv