08_Geoist之连续时间序列分析模块1

内容摘要：时间序列（TS）数据是最常见的一种地球物理测量数据表现形式。顾名思义，就是按照一定时间采样的测量数据。通常至少包括时间列和测量值两列数据。对于原始的时间序列，可能因为仪器故障、干扰、传感器漂移等问题，产生断记，曲线突然变换和长周期的单向漂移等问题。这时候，就需要对数据进行预处理。在此基础上，还可以进行一些分析和基于模型进行预测。今天我们结合Geoist的Snoopy模块，介绍一下时间序列分析的功能。

1、时间序列分析

时间序列包括按照时间间隔的递增 (以分钟、小时、天、周为单位等) 来分类的数据使用。时间序列信号非常广泛，每天温度的变化、股票的涨跌、地震波形记录等等。要对一个时间序列进行分析、建模和预测，首先要了解序列的性质。

一般对TS序列进行分析之前，要判断其稳定性和随机性。因为，平稳性（stationary）是很多分析模型和工具应用的基础。

基本分析逻辑是：

• 判断TS序列是否稳定；如果是可以选择下一步分析模型，如果否，则首先要使数据变的平稳；

• 平稳后的序列就可以对其进行建模，并根据实际观测数据来拟合模型参数；

• 有了模型后，进一步可以对信号的变化进行预测，以及对观测量进行正常和异常的判断。

对于实际的TS数据，由于时间序列数据的离散性，许多时间序列数据集都有一个周期性的以及/或着内置在数据中的趋势元素。时间序列建模的第一步是计算现有周期 (在平稳时间段内的周期性模式)以及/或着数据向上或向下移动的趋势。

平稳系列是指系列的平均值不再是一个有关于时间的函数。随着趋势数据的增加以及时间的推移，系列的平均值会随时间而增加或减少 (比如随着时间推移房价的稳步上升)。对于周期性数据，系列的平均值随周期波动 (比如每24小时中，温度的上升和降低)。

一般有两种方法可用于实现平稳：差分数据或线性回归。差分数据指的是，你计算了两个连续观测中的差异；而线性回归则是，你可以在模型中为了周期性组件采用二进制指示器变量。在我们决定使用哪种方法之前，首先你自己必须非常了解数据和其中的物理意义。

（1）ADF测试

ADF测试(Augmented Dickey-Fuller test)是一种检定时间序列平稳性的方法。本质上是一种平稳的单元根测试。一般适合于有明显的趋势项的时间序列检测问题。通常根据ADF检测结果的p-value可以判断时间序列信号的稳定性。

（2）ARIMA 建模

如果分析的TS序列中，有很多信息是具有自相关特点的，可以采用ARIMA 模型，该模型的全称为自回归集成移动平均值 (AutoregressiveIntegrated Moving Average或ARIMA)模型。ARIMA 模型包含了用于描述周期和趋势的参数 (例如，在一周中有几天使用了虚变量和差分)，还包含了自回归和/或移动平均数条件来处理数据中嵌入的自相关性。

在确定了最适合数据趋势和周期的模型之后，您还必须有足够的信息来生成较为准确的预测，这些模型的能力仍然是有限的，因为它们并没有考虑到在过去的一段时间内，兴趣变量本身的相关性。我们将这种相关性称为自相关，这在时间序列数据中是十分常见的。如果数据具有自相关性，就像我们所做的那样，那么可能会需要额外的建模来进一步改进基线预测。

2、TS序列的Python分析工具

在时间序列分析方面，我们主要围绕Pandas工具包进行开发，所以，我们先介绍一下Pandas。Pandas是Python生态中非常强大的一个数据分析包，它也是基于Numpy来开发的，非常适用于对结构化的数据集进行分析。

Pandas对TS序列分析方面有很多接口，但是绝不是仅限于分析TS序列。尤其是在对结构化数据文件的读取方面，Pandas有非常多好用和方便的函数。

接下来我们，就看看如何导入数据以及开始分析过程。老规矩，我们不多废话直接上代码。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import geoist.snoopy.tsa as tsa

# parameters for loading data
data_path = Path(tsa.__file__).parent
orig_file = Path(data_path,"data",'50002_1_2312.txt')
water_file = Path(data_path,"data",'water_level_res.txt')

# parameters for processing data
na_values = None

# load data
data=pd.read_csv(Path(orig_file),parse_dates=[[0,1]],header=None,delim_whitespace=True,index_col=[0],na_values=na_values)
data.index.name = 'time'
data.columns = ['origin_data']
ax=data.plot(figsize=(15,12),y=data.columns[0])
ax.set_xlabel('Date')
ax.set_ylabel('Value')

大家注意前面import的模块，在Geoist的Snoopy目录下，我们带了两个示例数据，orig_file和water_file，大家可以找到它们直接打开看看。

然后，就是读入数据，我们用到了pandas的read_csv函数，返回的对象实例叫data，我们type(data)一下，信息如下：

Out[46]: pandas.core.frame.DataFrame

DataFrame大家一定要记住，这就是Pandas里面大名鼎鼎的数据框类型，简单理解就是一个类似excel表的数据结构，与数组不同的是，每列之间的数据类型可以不一致。通过数据框类型的Plot接口直接可以画图，记住不用再通过matplotlib的pyplot来画图啦！上述代码运行结果如下：

图1原始的时间序列数据

图中最需要注意的就是2009，2010以及很多地方的突然跳跃，这些信号很多情况都与仪器的观测故障相关，因此，我们必须去掉这些问题。

3、尖峰、突跳修补

这时候就要Geoist登场啦，通过前面导入的tsa模块里面的despike_v2函数，对原始数据进行处理，这里面我们引入thresh_hold参数，这个参数值是判断是否为突跳的关键。

详细过程看下面代码：

# despike
thresh_hold = 200.0

data['despiked'],data['flag'] = tsa.despike_v2(data['origin_data'].interpolate(),th=thresh_hold)
ax=data.plot(figsize=(15,12),y=data.columns[:-1])
ax.set_xlabel('Date')
ax.set_ylabel('Value')
plt.grid()
plt.legend()
plt.title("The preliminary result by threshold={}".format(thresh_hold),loc='left')

图2是处理结果，橙色曲线是预处理过后的，这个样子是不是看起来更正常一些呢？

图2 去尖峰突跳后的曲线效果对比

4、ADF平稳性测试

TS信号没有数据错误后，就可以开始测试了。我们首先来一个ADF测试，同样还是tas模块中的功能，函数为adfuller。详细代码如下：

# test stationary
window_size = 50

res = tsa.adfuller(data['despiked'].values)
tsa.print_adf(res,'despiked data')
data['mean'] = data['despiked'].rolling(window=window_size).mean()
data['std'] = data['despiked'].rolling(window=window_size).std()
data.plot(figsize=(15,12),y=['despiked','mean','std'])

大家注意我们又新定义了一个window_size参数，这个参数是每次测量的TS序列窗长，给定不同的值，结果可能会有一定的变化。运行后，结果如图3所示，我们可以看到三条曲线。

图3 ADF平稳性测试

print_adf函数输出了量化的测试结果，我们可以注意一下p-value，因为，一会对TS信号处理后，还会继续进行测试，我们可以观察其变化。

Augmented Dickey-Fuller test for despiked data:
  adf: -0.8225755579593904
  p-value: 0.8124005630822038
  norder: 11
  number of points: 3903
  critical values:
    1% : -3.4320265580345004
    5% : -2.8622808115385583
    10% : -2.567164342737072

5、去趋势

上面测试结果中的p值较大，主要是与数据中的趋势项有关，那么，如果去掉趋势呢？我们可以用tsa中的detrend函数（这个detrend实质上就是曲线性），用法如下：

# detrend
data['detrend'] = tsa.detrend(data['despiked'])
data.plot(figsize=(15,12),y=['detrend','despiked'])

# test stationary again
res = tsa.adfuller(data['detrend'].values)
tsa.print_adf(res,'detrended data')
data['mean'] = data['detrend'].rolling(window=window_size).mean()
data['std'] = data['detrend'].rolling(window=window_size).std()
data.plot(figsize=(15,12),y=['detrend','mean','std'])

运行上述代码，结果如下：

图4 去趋势效果

再进行一次ADF测试，如图5所示。

图5 去趋势后的ADF测试结果

量化结果如下，对比一下p值是不是小多了，这说明信号相对平稳了不少。

Augmented Dickey-Fuller test for detrended data:
  adf: -2.8626974817975075
  p-value: 0.04986090928335389
  norder: 11
  number of points: 3903
  critical values:
    1% : -3.4320265580345004
    5% : -2.8622808115385583
    10% : -2.567164342737072

6、周期项分解

如果您对直接去线性这种方法不满意，感觉模型太简单。那些可以采用STL方法，该方法是一种把时间序列分解为趋势项(trend component)、季节项(seasonal component)和余项(remainder component)的过滤过程。

同样，Geoist的snoopy模块也提供了这个函数，在tsa模块的seasonal_decompose函数中，用法如下：

# seasonal decomposition
period = 365
na_values_output = np.nan
decomposition = tsa.seasonal_decompose(data['despiked'],freq=period,extrapolate_trend='freq')
fig=decomposition.plot()
fig.set_size_inches(15,8)
data['trend'] = decomposition.trend.fillna(na_values_output)
data['seasonal'] = decomposition.seasonal.fillna(na_values_output)
data['residual'] = decomposition.resid.fillna(na_values_output)

# test stationary on residual
res = tsa.adfuller(data['residual'].dropna().values)
tsa.print_adf(res,'residual data')
data['mean'] = data['residual'].rolling(window=window_size).mean()
data['std'] = data['residual'].rolling(window=window_size).std()
data.plot(figsize=(15,12),y=['residual','mean','std'])

做完STL分解后，我们又进行了一遍ADF测试。分解结果如图6所示。

图6 周期项分解

图7是ADF测试结果，大家可以对比一下均值和标准差（mean和std）的变化呦！

图7 再测试结果

量化结果如下，大家注意p值是不是又小了不少，这说明方法是有效的。

Augmented Dickey-Fuller test for residual data:
  adf: -5.1966236577048
  p-value: 8.932706805879938e-06
  norder: 13
  number of points: 3901
  critical values:
    1% : -3.432027418154056
    5% : -2.8622811914877873
    10% : -2.567164545006864

7、差分、移动平滑、指数平滑

TS信号的处理除了预处理，稳定性检测和去趋势外，还包括常用的差分、平滑等常规操作，还应该知道Pandas的dataframe数据类型所提供的这些接口：diff，rolling，ewm。
示例如下：

(1)差分

通常对时间序列求导，需要进行差分运算，这些常规的数据处理不需要自己再开发，直接利用dataframe的数据类型接口就行。

# first order difference
data['diff'] = data['origin_data'].diff()
ax = data.plot(figsize=(15,12),y=['origin_data','diff'])

# second order difference
data['diff2'] = data['diff'].diff()
ax = data.plot(figsize=(15,12),y=['origin_data','diff2'])

(2)滑动平均

滑动平均，最常用的一种低通滤波方法，通过rolling函数可以容易实现，方法如下：

# moving average
window_size = 50
center = False
data['ma'] = data['residual'].rolling(window=window_size,center=center,min_periods=1).mean()
data.plot(y=['residual','ma'])

图8 滑动平均滤波

该方法的窗口取得越长，噪声被去除的就越多，我们得到的信号就越平稳；

但同时，信号的有用部分丢失原有特性的可能性就越大，而我们希望发现的规律丢失的可能性就越大。

(3)指数平滑

除了滑动平均外，还有指数平滑法，该方法由布朗(Robert G..Brown)提出，他认为时间序列的态势具有稳定性或规则性，所以时间序列可被合理地顺势推延。

上式中，St为时间t的平滑值，yt为观测值，a为平滑常数。

该方法认为最近的过去态势，在某种程度上会持续到最近的未来，所以将较大的权数放在最近的资料。

指数平滑法是生产预测中常用的一种方法。也用于中短期经济发展趋势"预测")，所有预测方法中，指数平滑是用得最多的一种。

简单的全期平均法是对"时间数列"的过去数据一个不漏地全部加以同等利用；移动平均法则不考虑较远期的数据，并在"加权移动平均法"中给予近期资料更大的权重；

而指数平滑法则兼容了全期平均和移动平均所长，不舍弃过去的数据，但是仅给予逐渐减弱的影响程度，即随着数据的远离，赋予逐渐收敛为零的权数。

也就是说指数平滑法是在"移动平均法"基础上发展起来的一种 "时间序列分析预测法"，它是通过计算指数平滑值，配合一定的时间序列预测模型对现象的未来进行"预测"。

其原理是任一期的指数平滑值都是本期实际观察值与前一期指数平滑值的加权平均。

实际的用法如下：

# exponential moving average
factor = 0.3
data['ewm'] = data['residual'].ewm(alpha=factor).mean()
data.plot(y=['residual','ewm'])

上面代码factor取值范围为0-1，取值约接近于1，前面的信号对平滑结果影响越小。

图9 指数平滑效果

8、距平分析

距平是某一系列数值中的某一个数值与平均值的差，分正距平和负距平。距平分析常用于分析具有显著周期性的信号，如：气象观测数据。

距平值用来确定某个时段或时次的数据，相对于该数据的某个长期平均值。

举个例子：一个地区某天的平均气温是14度，该地区该天平均气温的30年平均值是12度，那么该地区该天的平均气温距平就是2度。2度的距平表明今天的平均气温相对于该地区该天平均气温的30年平均值偏高2度。但是人体所能感觉的真实气温是14度。

距平分析方法，由tsa模块的departure函数实现，用法如下：

# load data
dateparser = lambda x: pd.to_datetime(x,format='%Y%m')
water = pd.read_csv(water_file,header=None,parse_dates=True,index_col=0,delim_whitespace=True,date_parser=dateparser)
water[water == 99999] = np.nan
water = water.interpolate()
water.columns = ['origin','mean','departure']
water_origin = pd.DataFrame(water[water.columns[0]]).copy()

# call departure, and plot.
water_origin,_ = tsa.despike_v2(water_origin,th=200)
wate_departure = tsa.departure(water_origin)
ax = wate_departure.plot(figsize=(16,9))
ax.invert_yaxis()

距平分析结果如图10所示，图中绿色的departure曲线就是距平值，在2009-2010年的变化最大。

图10 矩平分析结果

9、变采样、月平均与滤波

Pandas的dataframe还支持重采样功能，resample函数，用法如下：

# upsample
water_daily = water_origin.resample('D').asfreq().interpolate()
water_daily.head(10)

# downsample
water_monthly = water_daily.resample('MS').asfreq().interpolate()
water_monthly.head(10)

# monthly mean
water_monthly = water_daily.resample('MS').mean().interpolate()
water_monthly.head(10)

说了这么半天，还都仅限于时间域的分析方法，没有说频率与的方法。没关系，tas模块也提供了滤波函数，方法如下：

# filter
# generate dataset
sample_rate = 30.0
n = np.arange(300)
orig_data = np.sin(0.1*np.pi*n)+2.0*np.cos(0.5*np.pi*n)+1.5*np.sin(0.8*np.pi*n)

# generate filter
order = 10
nyq = 0.5*sample_rate
lower_cut_rate = 7.0 / nyq
upper_cut_rate = 10.0 / nyq
sos = tsa.butter(10,lower_cut_rate,btype='low',output='sos')

# apply filter to data
filtered_data = tsa.sosfiltfilt(sos,orig_data)

# plot data
fig = plt.figure(figsize=(16,16))
ax = plt.subplot(211)
ax.plot(n/sample_rate,orig_data,label='orig_signal')
ax.plot(n/sample_rate,filtered_data,label='filtered_signal')
ax.set_xlabel('time(s)')
ax.set_ylabel('magnitude')
ax.legend()
plt.title('Effect of low pass filter (critical frequency:{}Hz)'.format(lower_cut_rate*nyq),loc='left')
plt.grid()

ax = plt.subplot(212)
w,h = tsa.sosfreqz(sos,worN=1000)
ax.plot(0.5*sample_rate*w/np.pi,np.abs(h))
ax.set_xlabel('frequency(Hz)')
ax.set_ylabel('response')
plt.title('Frequency response',loc='left')
plt.grid()

上边的信号是模拟产生的，效果如下：

图11 滤波和频谱

10、功率谱分析

既然谈到了频率分析方法，就离不开功率谱密度（PSD）。

因为物理学中，信号通常是波的形式表示，例如电磁波、随机振动或者声波。当波的功率频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）

记住welch、lombscargle和periodogram这三个函数即可，具体方法如下：

# psd
f_w,pxx_w = tsa.welch(orig_data,sample_rate,nperseg=256,scaling='spectrum')
f_p,pxx_p = tsa.periodogram(orig_data,sample_rate,scaling='spectrum')
f_l = np.linspace(0.1,14,3000)*np.pi*2.0
pxx_l = tsa.lombscargle(n/sample_rate,orig_data,f_l)

# plot result
fig = plt.figure(figsize=(15,9))
ax = fig.add_subplot(111)
ax.plot(f_w,pxx_w,label='welch')
ax.scatter(f_p,pxx_p,label='peridogram',c='g')
ax.plot(0.5*f_l/np.pi,np.sqrt(pxx_l*4.0/len(orig_data)),alpha=0.7,label='lombscargle')
ax.legend()

图12 功率谱分析

11、ARIMA建模

时间序列的预处理包括两个方面的检验，平稳性检验和白噪声检验。能够适用ARMA模型进行分析预测的时间序列必须满足的条件是平稳非白噪声序列。

关于ARIMA模型的介绍在文中开头就简单讲述了，下面我们直接看实现方法：

# ARIMA
p = 5
d = 0
q = 1
P = 0
D = 0
Q = 0
s = 0
model = tsa.SARIMAX(data['detrend'].dropna(),
                    order=(p,d,q),
                    seasonal_order=(P,D,Q,s),
                    enforce_stationarity=False)
rests = model.fit()

pred = rests.get_forecast(180)
pci = pred.conf_int()
fig, (ax0, ax1, ax2) = plt.subplots(nrows=3, ncols=1,figsize=(12,12))
ax0.plot(data['detrend'].dropna())
ax1.plot(data['detrend'].dropna().values)
pred.predicted_mean.plot(ax=ax1,label='forecast')
ax1.fill_between(pci.index,pci.iloc[:,0],pci.iloc[:,1],color='k',alpha=0.2,label='0.05 confidence interval')
ax1.legend()

start_day=365*10
pred = rests.get_prediction(start_day)
pci = pred.conf_int()
pred.predicted_mean.plot(ax=ax2,label='prediction')
ax2.fill_between(pci.index,pci.iloc[:,0],pci.iloc[:,1],color='k',alpha=0.2,label='0.05 confidence interval')
data['detrend'].dropna().iloc[start_day:].asfreq('5d').plot(style='o',ax=ax2)
ax2.legend()

其中，get_forecast为向外预测（Out-of-sample forecasts）；而get_prediction具有内外预测功能（In-sample prediction and out-of-sample forecasting），可以作为一种滤波器使用。

图13上图是去趋势后的结果，中间图最后为外推180天的结果和误差估计，下图为2018年4月的预测和实际观测结果对比。

图13 模型拟合与预测

一句话总结时间序列分析的根本还是去除信号中的“杂质”，而平稳过程信号，一般才具有可以预测的意义。信号中的随机干扰永远是要去除的部分，这也是模型能帮你做的事情。TS序列分析很复杂，模型和理论有很多，如果你对这方面感兴趣，可以进一步阅读状态空间模型和Kalman滤波方面的资料，这些都是实战中必不可少的基础知识。