0 前言
霍克斯过程在金融领域常被用作高频交易的算法,通过对毫秒级的订单簿的买卖行为进行建模,并通过结果构建交易策略。
这篇文章选用商品期货合力合约作为标的,面向霍克斯过程在分钟行情数据的应用探索, 时间范围跨越7个月。
由于交易所的存在休市,事件的到达之间存在一定的空白时间段,而这导致时间到达时间不再连续,这部分本文尚未处理。
一级目录 | 二级目录 | 三级目录 |
---|---|---|
霍克斯过程 | ||
实践的流程 | 方案构思 | |
探索性分析 | ||
数据建模 | ||
模型评估 | Q-Q 残差图 | |
日级分布对比 | ||
参数置信区间 | ||
月级参数推理 |
1 霍克斯过程
霍克斯过程是用于对自激励过程建模的一种计数过程,常常用于分析地震、流行病、社交网络、金融市场等的动态数据。
它描述在时间维度发生的一系列事件,这些事件的发生都会激励下一件事情发生的可能性,同时这个激励的效果会随着时间衰减。
它的核心概念是条件强度函数,它表示在给定过去发生的所有事件的情况下,在某个时间点附近发生事件的期望频率。
条件强度函数的数学表达式如下:
其中
是条件强度函数; 是第 个事件发生的时间; 是背景强度,表示在没有任何激励效应的情况下,事件发生的基础频率; 是触发核,表示每个过去的事件对当前时间点的激励强度。
触发核最常见的形式是指数衰减函数:
其中, 是触发核的幅度, 是衰减率。
指数衰减函数用于表示过去的事件对当前时间点的激励效果随着时间间隔的增加而出现指数的衰减。
2 实践的流程
2.1 方案构思
为了将霍克斯过程应用在商品期货分钟行情中,我们首先需要对事件进行定义,这里简单地将每一分钟间隔发生正收益率现象作为事件,简称为上涨事件。
在使用霍克斯过程对数据进行建模之前,我们通过探索性分析事先验证上涨事件的到达是否符合霍克斯过程的特性。
在建模后使用一系列检验方法来分析最后的模型的质量,以佐证模型拟合优度。
本文选取的分钟数据是螺纹钢主力合约 RB2405 的一分钟行情,采样范围是从 2023年5月16日 到 2023年12月29日,总共有 52277 条记录。
2.2 探索性分析
首先我们对数据进行描述性统计,这主要是来看各字段的异常和分布,可以看到数据中并没有特殊的离群值。(特别关注的是后续需要使用的收盘价字段)
df.describe()
open | high | low | close | amount | volume | |
---|---|---|---|---|---|---|
count | 52277 | 52277 | 52277 | 52277 | 5.23e+04 | 52277 |
mean | 3718.908545 | 3720.077644 | 3717.693345 | 3718.884289 | 5.08e-40 | 627.518469 |
std | 165.227356 | 165.493968 | 164.98057 | 165.253576 | 6.25e-40 | 1802.452899 |
min | 3321 | 3322 | 3320 | 3322 | 1.54e-44 | 0 |
25% | 3607 | 3608 | 3607 | 3608 | 4.59e-41 | 10 |
50% | 3685 | 3687 | 3684 | 3685 | 2.42e-40 | 56 |
75% | 3844 | 3845 | 3843 | 3844 | 6.20e-40 | 375 |
max | 4074 | 4076 | 4072 | 4075 | 2.07e-39 | 48057 |
根据每一分钟前后的收盘价计算收益率,并标记大于0为上涨事件。(简单认为在前一分钟末开仓,在后一分钟末进行平仓)
tmp = df.groupby(df.index.date).apply(lambda x: x['close'].pct_change() > 0).droplevel(level=0).to_frame()
tmp = tmp.rename(columns = {"close": "return"})
df = pd.merge(df, tmp, how='left', on='date')
df
date | open | high | low | close | amount | volume | return |
---|---|---|---|---|---|---|---|
2023-05-16 09:01:00 | 3565 | 3565 | 3552 | 3562 | 9.71e-43 | 9 | False |
2023-05-16 09:02:00 | 3562 | 3562 | 3562 | 3562 | 9.71e-43 | 0 | False |
... | ... | ... | ... | ... | ... | ... | ... |
对每一天的分钟行情数量进行统计,数量一直稳定在300条以上,认为缺失影响较小。(图中出现归零,是因为这些天并不在交易日期内)
t_data = df['close'].resample('d').count()
plt.plot(t_data.index, t_data.values)
为了方便分析,我们将分钟行情聚合成天级,并仅保留处于交易日行情。
p_data = df.resample('d').agg({'open': 'first', 'high': 'max', 'low': 'min', 'close': 'last', 'amount': 'sum', 'volume': 'sum', 'return': 'sum'})
p_data = p_data[p_data.index.map(lambda row: row in trade_date['trade_date'].tolist())]
p_data
date | open | high | low | close | amount | volume | return |
---|---|---|---|---|---|---|---|
2023-05-16 | 3565 | 3578 | 3516 | 3552 | 5.00e-40 | 2398 | 74 |
2023-05-17 | 3535 | 3599 | 3511 | 3532 | 8.61e-40 | 2405 | 78 |
... | ... | ... | ... | ... | ... | ... | ... |
在日级颗粒度的行情数据中,我们分别进行上涨事件到达数量和交易量的统计,并勾勒图形。(第一幅图在存在着趋势,这可能行情长期趋势导致的)
plt.plot(p_data.index, p_data['return'])
plt.plot(p_data.index, p_data['volume'])
2.3 数据建模
提取上涨事件的数据。(随便一个列都可以)
up_data = df[df['return']]
根据它的到达索引进行建模,并使用最大似然法获取模型内置参数估计值。
timeToQuake = up_data.index - up_data.index[0]
ts = np.array(timeToQuake.total_seconds() / 60 / 60 / 24)
obsPeriod = up_data.index[-1] - up_data.index[0]
T = obsPeriod.days * 252 // 365 * 4 * 60
𝛉_exp_mle = hawkes.exp_mle(ts, T)
print("Exp Hawkes MLE fit: ", 𝛉_exp_mle)
Exp Hawkes MLE fit: [1.21443093e-02 1.51569882e+02 1.55197132e+02]
2.4 模型评估
2.4.1 Q-Q 残差图
通过 Q-Q 残差图来检验残差数据是否近似于正态分布。(一般而言,模型拟合越好,残差会越接近正态分布)
tsShifted = hawkes.exp_hawkes_compensators(ts, 𝛉_exp)
iat = np.diff(np.insert(tsShifted, 0, 0))
qqplot(iat, dist=stats.expon, fit=False, line="45")
plt.show()
2.4.2 日级分布对比
将原始数据以及模型仿真的数据分别进行分布的比较。(如果模型拟合较好,那么它模拟生成的数据应会接近原始数据分布)
t_date = [timedelta(seconds=x) + up_data.index[0] for x in ts * 24 * 60 * 60]
t_df = pd.DataFrame({'date': t_date, 'count': np.ones(len(t_date))})
t_df.set_index('date', inplace=True)
a_data = t_df.resample('d').sum()
a_data = a_data[a_data.index.map(lambda row: row in trade_date['trade_date'].tolist())]
plt.plot(a_data.index, a_data['count'])
tsSim = hawkes.exp_simulate_by_thinning(𝛉_exp, len(ts))
t_date = [timedelta(seconds=x) + up_data.index[0] for x in tsSim * 24 * 60 * 60]
t_df = pd.DataFrame({'date': t_date, 'count': np.ones(len(t_date))})
t_df.set_index('date', inplace=True)
t_df.resample('d').sum().plot()
2.4.3 参数置信区间
使用上述求解的参数多次进行仿真采样,并使用这些采样的数据进行重新建模。(重复迭代一百次,会得到一百组参数,使用它们的百分位可以看到参数的波动区间)
θ_bootstrap = []
for seed in range(100):
rnd.seed(seed)
tsSim = hawkes.exp_simulate_by_thinning(𝛉_exp, T)
θ_mle = hawkes.exp_mle(tsSim, T)
θ_bootstrap.append(θ_mle)
θ_bootstrap = np.vstack(θ_bootstrap)
print(np.quantile(θ_bootstrap, [0.025, 0.975], axis=0))
[[1.09225173e-02 1.31199006e+02 1.46149548e+02] [1.31185051e-02 1.57786689e+02 1.64188471e+02]]
2.4.4 月级参数推理
将上涨事件拆分成月级,每一个月进行一次参数估计,以评估参数随时间的平稳性。(有些类似模型外推的泛化性检验)
months = sorted(set(up_data.index.month))
numMonths = len(months)
θ_monthly = []
for i, month in enumerate(months):
returnMonth = up_data[up_data.index.month == month]
timeToReturn = returnMonth.index - returnMonth.index[0]
ts_i = timeToReturn.total_seconds() / 60 / 60 / 24
obsPeriod = returnMonth.index[-1] - returnMonth.index[0]
T_i = obsPeriod.days * 252 // 365 * 4 * 60
θ_monthly.append(hawkes.exp_mle(np.array(ts_i), T_i))
θ_monthly = np.vstack(θ_monthly)
plt.plot(θ_monthly)
plt.show()
3 参考文章
《霍克斯Hawkes过程》:https://www.jianshu.com/p/5c40fbd8eaf1 《Hawkes Process:超详细带示例的讲解》:https://blog.csdn.net/weixin_44835327/article/details/130022377 《霍克斯过程介绍以及一维霍克斯过程求解》:https://zhuanlan.zhihu.com/p/566644881 《霍克斯过程Hawkes Processes分析比特币交易数据订单到达自激过程时间序列》:https://mp.weixin.qq.com/s/MhUVHyeCQy2wt7IvR7IqFg 《The Elements of Hawkes Processes》:https://link.springer.com/book/10.1007/978-3-030-84639-8 《Bitcoin Trade Arrival as Self-Exciting Process》:http://jheusser.github.io/2013/09/08/hawkes.html