Python 之 时间序列测试(TimeSeriesTests)

时间序列是一个随机变量的序列。

时间序列是描述一个在时间上发生的过程的数值序列，在连续的时间点上测量，通常是以相等的时间间隔。

时间结构给观察结果增加了一定的秩序。这个顺序意味着，关于现有观测结果一致性的重要假设需要单独处理。

例如，检查一行的静止性。

在时间序列中，通过加入趋势性、季节性等与时间相关的结构，可以很容易地打破对时间序列分布参数的恒定性（静止性）的假设。

如果我们想保证一个固定模型，我们假设我们的数据是一个固定过程的实现。因此，我们分析的第一步应该是检查数据中是否有趋势或季节性影响的迹象，如果有，就将其删除。

检查静止性:

图表
统计摘要
静止统计测试

---

目录

1.导入所需模块和数据

2.静态/固定时间序列

2.1视觉检测​

2.2时间序列构成部分的分配

3.统计检验

3.1固定性

3.2规范性​

3.2同源性

4.不稳定行

4.1图形检测​

4.2时间序列构成部分的分配

4.3统计检验​

5.下载数据dbnomics​

---

# settings
alpha = 0.05

1.导入所需模块和数据

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.stats.stattools import jarque_bera
from scipy import stats
from statsmodels.iolib.table import SimpleTable

import itertools

series = pd.read_csv('data/daily-total-female-births.txt', sep=',', header=0, index_col=0)
series.index = pd.to_datetime(series.index)  # for series decomposition

2.静态/固定时间序列

2.1视觉检测

series.plot()
plt.show()

resampled_series = series.resample('W', how='mean')
resampled_series.plot()
plt.show()

大家可以看到，新图有一个亮丽的展示（9-10月），没有明显的趋势。

series.hist()
plt.show()

我们清楚地看到高斯分布的钟形，可能有较长的右尾。

result_hist = series.describe()
result_hist

从特征图和直方图中可以看出，该序列的特征比较均匀，变异比较小，由变异系数可以看出：𝑉=𝜎𝑥¯，其中𝜎为标准差，𝑥¯为样本的算术平均值。在我们的情况下，它是平等的。

print('V = %f' % (result_hist.loc['std', 'Births'] / result_hist.loc['mean', 'Births']))

split = int(series.shape[0] / 2)
series1, series2 = series.values[:split], series.values[split:]
mean1, mean2 = series1.mean(), series2.mean()
var1, var2 = series1.var(), series2.var()
print('mean1=%f, mean2=%f' % (mean1, mean2))
print('variance1=%f, variance2=%f' % (var1, var2))

一个固定/静态过程的例子：白噪声。

2.2时间序列构成部分的分配

decomposition_result =seasonal_decompose(series)

decomposition_result.plot()
plt.show()

series.plot.box()
plt.show()

3.统计检验

3.1固定性

from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.stattools import kpss

adfuller_result = adfuller(series.squeeze())
print('ADF Statistic: %f' % adfuller_result[0])
print('p-value: %f' % adfuller_result[1])
print('lags used: %f' % adfuller_result[2])
print('Critical Values:')
for key, value in adfuller_result[4].items():
    print('\t%s: %.3f' % (key, value))

结果解读

作为输出的一部分，我们得到一个搜索表，帮助我们确定ADF统计。我们可以看到，我们的值-4.8小于-3.449，对应的p值为1%。

这说明我们可以用低于1%的值（即结果是统计学上的巧合的概率很低）拒绝零假设。

拒绝零假设意味着过程没有单一的根，进而意味着时间序列是静止的或不具有时间依赖性结构。

kpsstest = kpss(series.values.squeeze(), regression='c')
print('KPSS Statistic: %f' % kpsstest[0])
print('p-value: %f' % kpsstest[1])
print('lags used: %f' % kpsstest[2])
print('Critical Values:')
for key, value in kpsstest[3].items():
    print('\t%s: %.3f' % (key, value))

根据得到的KPSS检验统计量值，拒绝时间序列非平稳性的零假设。

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

plot_acf(series.values.squeeze(), lags=25)
plt.show()

plot_pacf(series.values.squeeze(), lags=25)
plt.show()

3.2规范性

k2, p = stats.normaltest(series)
print("p = {}".format(round(p[0], 5)))
if p < alpha:  # null hypothesis: series comes from a normal distribution
    print("Нулевая гипотеза может быть отвергнута.")
else:
    print("Нулевая гипотеза не может быть отвергнута, выборка - из нормального распределения.")

3.2同源性

当正态分布有显著偏差时，Levene检验是Bartlett检验（scipy.stats.bartlett）的替代方法。

series_parts = np.array_split(series.values.squeeze(), 4)
index_combinations = itertools.product(list(range(4)), list(range(4)))
count_homosced = 0
count_heterosced = 0
for idx1, idx2 in index_combinations:
    statistic, p = stats.levene(series_parts[idx1], series_parts[idx2])
    if p < alpha:  # null hypothesis: all input samples are from populations with equal variances
        count_heterosced += 1
    else:
        count_homosced += 1

if count_homosced > count_heterosced:
    print("Нулевая гипотеза не может быть отвергнута - присутствует гомоскедастичность.")
else:
    print("Нулевая гипотеза может быть отвергнута - присутствует гетероскедастичность.")

4.不稳定行

air_passengers_series = pd.read_csv('data/airline-passengers.txt', sep=',', header=0, index_col=0)
air_passengers_series.index = pd.to_datetime(air_passengers_series.index)  # for series decomposition

4.1图形检测

air_passengers_series.plot()
plt.show()

air_passengers_series.hist()
plt.show()

result_hist = series.describe()
result_hist

print('V = %f' % (result_hist.loc['std', 'Births'] / result_hist.loc['mean', 'Births']))

air_passengers_log_series = np.log(air_passengers_series.values)
plt.plot(air_passengers_log_series)
plt.show()

plt.hist(air_passengers_log_series)
plt.show()

split = int(air_passengers_series.shape[0] / 2)
series1, series2 = air_passengers_series.values[:split], air_passengers_series.values[split:]
mean1, mean2 = series1.mean(), series2.mean()
var1, var2 = series1.var(), series2.var()
print('mean1=%f, mean2=%f' % (mean1, mean2))
print('variance1=%f, variance2=%f' % (var1, var2))

4.2时间序列构成部分的分配

decomposition_result = seasonal_decompose(air_passengers_series)
decomposition_result.plot()
plt.show()

我们看到了明显的趋势性和季节性。

air_passengers_series.plot.box()
plt.show()

4.3统计检验

adfuller_result = adfuller(air_passengers_series.squeeze())
print('ADF Statistic: %f' % adfuller_result[0])
print('p-value: %f' % adfuller_result[1])
print('lags used: %f' % adfuller_result[2])
print('Critical Values:')
for key, value in adfuller_result[4].items():
    print('\t%s: %.3f' % (key, value))

kpsstest = kpss(air_passengers_series.values.squeeze(), regression='c')
print('KPSS Statistic: %f' % kpsstest[0])
print('p-value: %f' % kpsstest[1])
print('lags used: %f' % kpsstest[2])
print('Critical Values:')
for key, value in kpsstest[3].items():
    print('\t%s: %.3f' % (key, value))

不能拒绝单根的零假设（因此是非稳态性）。

plot_acf(air_passengers_series.values.squeeze(), lags=25)
plt.show()

plot_pacf(air_passengers_series.values.squeeze(), lags=25)
plt.show()

非稳态过程的例子：具有确定性趋势的过程，随机游走。

摆脱不稳定的状态：

以系列的差异。
以该系列的季节性差异为例
转换（对数化）

5.下载数据dbnomics

import pandas as pd
from dbnomics import fetch_series, fetch_series_by_api_link

df1 = fetch_series('AMECO/ZUTN/EA19.1.0.0.0.ZUTN')

df1.head()

以下展示部分截图

df1 = df1.dropna()
df1.head()

df2 = fetch_series_by_api_link("https://api.db.nomics.world/v22/series/WB/DB?observations=1&dimensions=%7B%22country%22%3A%5B%22FR%22%2C%22IT%22%2C%22ES%22%5D%2C%22indicator%22%3A%5B%22IC.REG.COST.PC.FE.ZS.DRFN%22%5D%7D")
df2.head()