使用希尔伯特-黄变换（HHT）进行时间序列分析

量化投资与机器学习微信公众号

发布于 2019-02-26 16:26:39

3.1K0

发布于 2019-02-26 16:26:39

文章被收录于专栏：量化投资与机器学习

本期作者：Bradbeer

本期编辑：阿威 | 公众号翻译部成员

将非平稳时间序列用经验模态分解（EMD）转为固有特征方程式并且捕获其趋势。可以尝试使用HHT，当然这只是其中的一种方法，并没有像其他方法一样存在数学证明等。

数据准备

为了方便起见，我们选取了富时100指数（FTSE100）过去10年的收盘价作为金融时间序列，只是作为我们研究，大家可以用其他指数。之后，我们会选取希尔伯特谱来分析固有特征方程式来提取即时数据信息。

富时100指数数据的提取式这样的，加载到dataframe里：

时间,   开盘价,   收盘价,   最高价,    最低价,  成交量
02-Jan-2009,4434.20,4561.80,4561.80,4430.00,407295392
05-Jan-2009,4561.80,4579.60,4618.10,4520.80,836675968
…
21-Dec-2018,6711.93,6721.17,6733.00,6653.65,1636792576
24-Dec-2018,6721.17,6685.99,6721.17,6661.04,173149664

import os
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
data_directory = 'some path...'

data_file = os.path.join(data_directory, 'ftse100_2009_2018.csv' )

ftseinfo = pd.read_csv(data_file)

测试数据

# convert strings to dates
ftseinfo['Date'] = pd.to_datetime(ftseinfo['Date'], format="%d-%b-%Y" )

# convenience variable for plots
date_axis = ftseinfo['Date']

# The date entries are not evenly spaced, so calculate number of days from first date for each
first_date = date_axis[0]
deltas = ftseinfo['Date'] - first_date

def getdays(delta):
    return delta.days

ftseinfo['deltadays'] = deltas.apply(getdays)

# convenience variable for days
days = ftseinfo['deltadays']
# plot closing price
fig = plt.figure(figsize=(12,6))
sp1 = fig.add_subplot(111)
plt.xlabel('date')
plt.ylabel('price')

sp1.plot_date( x=date_axis, y=ftseinfo['Close Price'] )

plt.legend(loc='center right')
plt.show()

当股价是非平稳状态并且有随机走势基础时，以上图形可以看作有线性的趋势结构，让我们加以测试。

平稳性

一组具有平稳性的时间序列不以某些时间点为依赖，并且没有趋势性或季节性。股价按理来说应该是平稳性的，可是很显然，有些走势并非如此。它可能具备一些趋势。当然我们不能肯定未来的价格依赖于过去的价格，趋势只是过去时间点已知发生的事实。所以，给出一段股价走势，这些数据可能呈现出非平稳定性。为了更好地分析，我们要去除这一非平稳状态。通常，我们通过观察数据差异（例如价格变动）而不是绝对价格。

Augmented Dickey-Fuller Test（用于测试稳态）：

# Dickey-Fuller test
from statsmodels.tsa.stattools import adfuller

def adf_test(timeseries):
    dftest = adfuller(timeseries, autolag='AIC')
    print('ADF Statistic: %f' % dftest[0])
    print('p-value: %f' % dftest[1])
    print('Critical Values:')
    for key, value in dftest[4].items():
        print('\t%s: %.3f' % (key, value))

结果显示于下方：

adf_test( ftseinfo['Close Price'] )
ADF Statistic: -2.129286
p-value: 0.232903
Critical Values:
 1%: -3.433
 5%: -2.863
 10%: -2.567

我们用一个0.05的p-value作为阀值，很显然这里的p-value大于0.05，所以我们认为这组时间序列是非平稳的，需要改进。

分解与残差（趋势）

在下面的例子里，用pyhtt包来处理HHT分解：

from pyhht.utils import inst_freq
from pyhht import EMD

def decompose( x, y ):
    decomposer = EMD( y )
    imfs = decomposer.decompose()
    return imfs
def emd_analysis( days, data_values ):

    # decompose time series into intrinsic mode functions
    imfs = decompose(days, data_values)

    # extract the residue (overall trend)
    imf_residue = imfs[len(imfs)-1]

    return (imfs, imf_residue )
# Do the decomposition on the closing price time series
( imfs, imf_residue ) = emd_analysis( days, ftseinfo['Close Price'])
# plot the residue to see if a clear trend in there
fig = plt.figure(figsize=(12,6))
sp1 = fig.add_subplot(111)
plt.xlabel('date')
plt.ylabel('residue')
sp1.plot_date( x=date_axis, y=imf_residue, color='red' )
plt.show()

残差显示了非常明显的趋势性，从分解结果可以看出一条十分明显的趋势线，因此需要去除。之后我们在用ADF（Augmented Dickey-Fuller Test）来测试是否稳态。

趋势去除第一步

ftseinfo['trend_adjusted_1'] = ftseinfo['Close Price'] - imf_residue
# Let's look at adjusted prices
fig = plt.figure(figsize=(12,6))
sp1 = fig.add_subplot(111)
plt.xlabel('date')
plt.ylabel('trend adjusted price')

sp1.plot_date( x=date_axis, y=ftseinfo['trend_adjusted_1'], color='green' )
plt.axhline(0, color='black')

plt.legend(loc='lower right')
plt.show()

再用ADF测试：

# Dickey-Fuller test for stationary or not
adf_test( ftseinfo['trend_adjusted_1'] )
ADF Statistic: -3.329548
p-value: 0.013608
Critical Values:
 1%: -3.433
 5%: -2.863
 10%: -2.567

发现p-value小于0.05，表明已经是平稳时间序列了。

分解步骤2

# Do the decomposition on the price movement series
( imfs, imf_residue ) = emd_analysis( days, ftseinfo['trend_adjusted_1'])

现在测试剩余残差：

# Check residue to ensure no significant trend information remains
fig = plt.figure(figsize=(12,6))
sp1 = fig.add_subplot(111)
plt.xlabel('date')
plt.ylabel('residue')
sp1.plot_date( x=date_axis, y=imf_residue, color='red' )
plt.axhline(0, color='black')
plt.show()

看起来也是需要去除的。

趋势去除第二步

ftseinfo['trend_adjusted_2'] = ftseinfo['trend_adjusted_1'] - imf_residue
# Dickey-Fuller test for stationary or not
adf_test( ftseinfo['trend_adjusted_2'] )
ADF Statistic: -3.343714
p-value: 0.013034
Critical Values:
 1%: -3.433
 5%: -2.863
 10%: -2.567

我们再分解一次。

( imfs, imf_residue ) = emd_analysis( days, ftseinfo['trend_adjusted_2'])
# Check residue to ensure no significant trend information remains
fig = plt.figure(figsize=(12,6))
sp1 = fig.add_subplot(111)
plt.xlabel('date')
plt.ylabel('residue')
sp1.plot_date( x=date_axis, y=imf_residue, color='red' )
plt.axhline(0, color='black')
plt.show()

从图形看还是剩余了一些趋势，但是在十年中的增长很少，可以认为趋势基本不存在了。

复权价格EMD

从以上看来，HHT/EMD方法找到了两个趋势，我们都已经将他们去除了，剩余的数据已经没有了趋势。因此我们可以开始测试固有特征方程了。

# Plot each IMF in turn; the highest frequency IMFS are first and lowest are last
for i in range(0,len(imfs)-1):
    fig = plt.figure(figsize=(12,6))
    sp1 = fig.add_subplot(111)
    plt.xlabel('date')
    plt.ylabel('imf %d'%i)
    sp1.plot_date( x=date_axis, y=imfs[i], color='green' )
    plt.axhline(0, color='black')
    plt.show()

随着分解过程的进行，先去掉最高频率的函数，然后识别出频率越来越低的函数。