算法交易实战日记（九）—— 再次探索HMM在价格状态预测上的应用

大嘎好！时隔两月我又回来了！
好久没写文章了，刚好最近做隐马尔可夫的二次探索有了新的发现和理解，和大家分享一下。本文含有大量分析源码，感兴趣的朋友也可以自己用其它模型尝试分析～
$ ?/ [4 s; c, m1 R; i
( Z" u( @3 T6 ?
+ w+ [% l! N5 O; e4 n6 W8 w$ f大家都知道西蒙斯就是靠HMM做择时一炮打红的，隐马尔可夫模型作为一个经典的概率模型可解释性很强，可探索性也很强，它不需要复杂的特征工程，直接从价格波动本身出发就可以做分析。
7 s* }& D/ F  m( M有关隐马尔可夫模型的原理本文不详述，简单来说就是你能看到的每一个状态（显性状态）其实背后都有一个对应的隐性状态，而隐形状态之间是依据马尔可夫过程进行状态转移的，比如显性状态是今天你喜欢的姑娘没出门，隐性状态是因为今天下雨了，那么你通过天气预告知道明天不下雨的概率（隐性状态转移概率），那你就可以大概估摸着到底哪天去约你喜欢的姑娘出去吃饭了。
HMM的假设是比较理想的，也就是隐性状态是可枚举而且独立的，在金融序列中，人们已经人为地确定过一些价格状态，比如拉盘/砸盘/横盘等，那么在这个基础上，我们是否可以假定知道价格的所有状态，然后对应的去寻找隐性状态的转移概率呢？
<hr/>1. 准备工作
8 ]+ ?5 s% m" D- o* J那么首先导入需要的库和数据提取，我使用的数据集为ETH/USDT永续合约5min级bar级数据，数据的时间范围为2020-09-26 08:00:00 到 2021-08-11 20:25:00，共计数据条数为92022条。
此外，本文中使用的隐马尔可夫模型开源库为hmmlearn。
import numpy as np
import pandas as pd
. {+ L7 Q; j4 A- timport matplotlib.pyplot as plt
import warnings
7 w; }- y' O6 g9 F' d/ mwarnings.filterwarnings(&#39;ignore&#39;)
from tqdm import tqdm
/ n3 P, M9 S1 q' k9 G# ]/ k% \6 G, u8 aimport random
import datetime
# I  ]) `/ f! }- o$ Ifrom hmmlearn.hmm import GaussianHMM,GMMHMM,MultinomialHMM
from sklearn.model_selection import train_test_split
import itertools
/ y: t. P+ J  H* E
' M- v0 B" Q* y: w% d: l1 p6 o$ Tdf = pd.read_csv(&#39;5m.csv&#39;)2. 特征提取
0 A+ V9 y. U. b$ }虽然不需要做复杂的特征工程，但是简单的特征工程还是要做滴！毕竟要做价格状态的判断，单纯依据价格本身是很不合理的，而bar级数据中包含了OHLC四维特征，从价格本身和波动率来判断的话，可以做以下三种不同特征：
bar级价格变化率: (close-open)/open
bar级上涨随机波动: (high-open)/open
bar级下跌随机波动: (low-open)/open
那么就可以定义一个提取特征的函数了：  
def _extract_features(data):
    open_price = np.array(data[&#39;开盘价&#39;])
2 A4 n5 r4 t1 B7 G    close_price = np.array(data[&#39;收盘价&#39;])
0 y' p& d. u/ N6 D    high_price = np.array(data[&#39;最高价&#39;])
    low_price = np.array(data[&#39;最低价&#39;])

0 {# U/ w. \8 O# d2 U( o* r% {    # 计算收盘价、高价和低价的分数变化
    frac_change = (close_price - open_price) / open_price
    frac_high = (high_price - open_price) / open_price
    frac_low = (low_price - open_price) / open_price
" a) ], ]5 z' [/ R! w" r  }) Z
    return np.column_stack((frac_change, frac_high, frac_low)),frac_change, frac_high, frac_low3. 模型建立与状态解释
# I6 S' x  \3 v我们假设每个隐性状态对应的显性状态都是正态分布的，再初步估计市场一共有5个状态：  
$ c, a' m$ \' {' x9 J% W7 @_states = 5
. C1 x  ]% L! Z% ktrain_data, test_data = train_test_split(df, test_size=0.1, shuffle=False)
features, frac_change, frac_high, frac_low = _extract_features(train_data)
model = GaussianHMM(n_components=_states)
+ h( M, U: K7 y
- K& M7 D5 E4 z+ g3 S" A* ?model.fit(features)

## GaussianHMM
print(&#34;隐藏状态的个数:&#34;, model.n_components)
7 j( w- G0 k/ b6 M8 |; Yprint(&#34;状态转移矩阵&#34;)
+ R& C5 L# a* m/ [6 bprint(model.transmat_)
0 ]% s# x, I# L2 s0 [: d' \) Eprint(&#34;均值矩阵&#34;)
print(model.means_)
print(&#34;方差矩阵&#34;)
print(model.covars_)
. U, x$ ]. a) |3 N4 B+ F' w0 ~6 s# v
4 R0 D1 H8 q" Ntrans_mat = model.transmat_
! E3 w1 D. i! xmeans_mat = model.means_
covar_mat = model.covars_其中trans_mat为隐性状态转移概率矩阵，means_mat为每个隐性状态对应均值矩阵，cover_mat为每个隐性状态对应方差矩阵。
结果如下：  
隐藏状态的个数: 5
5 I: P! \* R) ?
8 K) V) C  d, i$ r* b状态转移矩阵
[[4.92841919e-01 3.94296981e-02 4.88723685e-02 4.18168766e-01
4 l" ~9 s( m. G7 X  6.87248292e-04]
: K+ Q0 ^+ H( Z [1.22888498e-01 3.54667807e-01 4.04748939e-01 7.50435445e-02
  4.26512115e-02]
[1.20944857e-01 3.16391469e-01 3.37760871e-01 1.97331090e-01
  2.75717137e-02]
5 k# y4 o1 L2 M! g, Q# b) v( Y  D8 ` [5.02057408e-01 6.09745059e-02 5.49501647e-02 3.80524622e-01
0 m" X& K. T- G6 ?9 O& m- T2 O  1.49329942e-03]
& r: E  o9 g3 S! E& ] [5.46560440e-07 1.66771635e-01 2.20357886e-01 1.45143404e-06
  6.12868481e-01]]
; K! s) e3 v& N* y5 x( l& u7 y- K! V
均值矩阵
4 w  }3 O$ K0 \. E& M2 F[[ 0.00100483  0.00200838 -0.00069315]
3 G+ F+ V  v- v5 n [-0.0042892   0.00139075 -0.00662071]
0 `7 c  J, b( Z1 f+ I$ q [ 0.00430167  0.00626541 -0.00129704]
[-0.00127088  0.00068099 -0.0023776 ]
[-0.0002794   0.01021583 -0.01188547]]

方差矩阵
% y$ h4 j6 S6 j/ m8 F$ y9 h' y. ^[[[1.78831746e-06 0.00000000e+00 0.00000000e+00]
. {% O2 E, K+ K, O/ g- z  [0.00000000e+00 1.79026813e-06 0.00000000e+00]
  [0.00000000e+00 0.00000000e+00 7.92414902e-07]]
; s$ a' C; y2 E; W: E
- m: x- _; F! s+ s, a! f [[9.22402278e-06 0.00000000e+00 0.00000000e+00]
( g% X$ Z  {$ f3 u& i/ \  [0.00000000e+00 3.06631445e-06 0.00000000e+00]
  [0.00000000e+00 0.00000000e+00 8.07040990e-06]]
' S# M+ |9 U7 Q7 i& Z: g( L  o
[[8.66631732e-06 0.00000000e+00 0.00000000e+00]
  [0.00000000e+00 6.85396774e-06 0.00000000e+00]
  [0.00000000e+00 0.00000000e+00 3.04307130e-06]]
! o; e. r, [6 n& a1 ]
) G+ W# h1 v. S' m- Y. } [[1.90478192e-06 0.00000000e+00 0.00000000e+00]
* R( K# K8 z7 k* p" A  [0.00000000e+00 8.76364129e-07 0.00000000e+00]
& f" E% |3 L/ ~& s  [0.00000000e+00 0.00000000e+00 1.96545207e-06]]

[[2.98827765e-04 0.00000000e+00 0.00000000e+00]
  [0.00000000e+00 2.25462090e-04 0.00000000e+00]
  [0.00000000e+00 0.00000000e+00 2.71001164e-04]]]通过均值矩阵可以大概判断分别是怎么样的五个状态：
; ]9 {' s/ U& k2 S6 e

; A0 n& F/ x  b. H  L通过均值矩阵可以大致将以上几种状态分别归类为（对应状态0-4）：震荡向上、跌得飞起、涨得飞起、震荡向下、上下乱抖。
为了证明上述分类准确性，查看一下训练集分类的价格波动分布：  
plt.figure(figsize=(15,10))
7 I0 C1 b2 d) p! H* |: ^! Afor i in tqdm(range(_states)):
    idx = (decode == i)
    plt.plot_date(train_data[&#39;开盘时间_&#39;].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;))[idx],frac_change[idx],&#39;.&#39;,label=&#39;%dth hidden state&#39;%i,lw=0.5)
    plt.legend()
    plt.grid(1)

（中间颜色看起来好像汉堡王啊 =_= 饿了）


& ]8 U4 L. v. q2 q, V7 _通过波动分布感觉分类较为准确，下面观察一下各个状态的转移概率矩阵：
5 f8 z7 j; c  M* F( I

显性分布也较为符合认知。  
4. 显性状态分布
知道了隐性状态的转移概率后，下面可以分别观察每个隐性状态对应的显性状态（每个bar的价格变化比例）的分布：  
  d/ x8 f4 W$ g7 Z8 ]
) F7 T# h) v$ o1 z4 d7 B( l6 m

# V6 p2 f! q0 l) V$ b" }. i) c涨得飞起
6 u' ~& p. P( v0 u+ d/ O
: }) e8 c1 M, _3 [" ^

+ D& z2 V) _' D. J) [震荡向上
0 v6 {4 D; b, ]' [. }. e: ~2 m
, Q6 Q- |& b2 t8 r
& O+ H3 }& `/ |
震荡向下
( j9 ^/ `; \: w! D2 C9 c- L
8 S# k5 x! Y! Y! Q( u
' W: P( D3 i% i) L/ s2 `
. n+ k9 k4 C$ h: A( Q跌得飞起
+ i$ a  D4 X: P2 ^& x
0 L. A. V  r* ]- B/ L( r
$ Q3 [& p' Y' T( b3 e3 `
上下乱抖
/ |+ n. f- F+ q3 Q, C
观察到涨得飞起和跌得飞起都是高峰度的正态分布，而震荡向上和震荡向下都是明显有偏的偏态分布，上下乱抖则有两个尖峰，且非常薄尾。
3 _# ]' w; q: Z  K( \通过以上几个显性状态的分布可以更精确的进行抽样。
5. 蒙特卡洛模拟
最后使用测试集来进行隐性状态转移模拟和对比，我们对接下来的100个时间戳进行蒙特卡洛模拟，以训练集的最后一个时间戳的状态为起始状态，模拟10000次，再从训练集的显性分布中进行对应抽样（这里我的抽样逻辑不太严谨，是假设的正态分布进行的抽样，实际上可以进行分布拟合再进行抽样或者直接从原数据中进行抽样都会更好一些），观察上下三个sigma之间的模拟结果：
/ M  [' p. g/ D0 O# }0 Qdef simulation(init_state,_states,num_simulations,periods,trans_mat,means_mat,covar_mat):
    states = [i for i in range(_states)]
8 f) E* d2 H. r    mean_arr = [means_mat[j][0] for j in range(_states)]
    sd_arr = [np.sqrt(covar_mat[0][0]) for i in range(_states)]
+ F* U. n1 b3 I& ]/ M+ w: }6 a* m' u5 Q   
    midline_trend = []
    up_one_sigma_trend = []
    up_two_sigma_trend = []
' t: q9 B  Y" A/ w    up_three_sigma_trend = []
    down_one_sigma_trend = []
' u! \% q: T$ U    down_two_sigma_trend = []
8 l/ w( }  d" i5 E. d8 U    down_three_sigma_trend = []
0 M' X- s4 d0 _$ L( k) j2 V    sims = []
* B/ f% S1 p! g8 z+ E' n7 R3 t    for j in tqdm(range(num_simulations)):
/ N0 ]* N" s8 R$ @' E& \/ u( A        sim = [0]*(periods+1)
1 B- K4 H3 q8 N+ H9 `0 l        sim[0] = init_state
1 b$ q% x$ q! M        #sim[0] = model.decode(features)[1][-1]
        trans_df = pd.DataFrame(trans_mat)
, R" X. t+ }8 X+ s2 f        for i in np.arange(1,periods+1):
) u  n$ D' b/ ]            sim = np.random.choice(states, 1, p = trans_df.loc[sim[i-1],:])[0]
        sim = sim[1:]
        sims.append(sim)
+ ]) J+ y- }' ~, C5 p# W
& y: ^! E" ^& t* ^        midline_trend.append([mean_arr[sim] for i in range(periods)])
- o) `( Y8 Y2 }        up_one_sigma_trend.append([mean_arr[sim]+1*sd_arr[sim] for i in range(periods)])
- d" J0 E1 E: V2 [, I8 ?        up_two_sigma_trend.append([mean_arr[sim]+2*sd_arr[sim] for i in range(periods)])
5 Z$ O. j; Z, E7 R; g        up_three_sigma_trend.append([mean_arr[sim]+3*sd_arr[sim] for i in range(periods)])
        down_one_sigma_trend.append([mean_arr[sim]-1*sd_arr[sim] for i in range(periods)])
        down_two_sigma_trend.append([mean_arr[sim]-2*sd_arr[sim] for i in range(periods)])
        down_three_sigma_trend.append([mean_arr[sim]-3*sd_arr[sim] for i in range(periods)])
! H  G5 ]7 [. Y+ P, B" c. ~        
    sims = pd.DataFrame(sims).mode()
    midline_trend = pd.DataFrame(midline_trend).mean()
    up_one_sigma_trend = pd.DataFrame(up_one_sigma_trend).mean()
/ _, w: o$ ?8 }7 X2 C; s5 P- X    up_two_sigma_trend = pd.DataFrame(up_two_sigma_trend).mean()
    up_three_sigma_trend = pd.DataFrame(up_three_sigma_trend).mean()
& B0 z5 @& }7 o    down_one_sigma_trend = pd.DataFrame(down_one_sigma_trend).mean()
+ X! _0 O! \! `( _" P    down_two_sigma_trend = pd.DataFrame(down_two_sigma_trend).mean()
    down_three_sigma_trend = pd.DataFrame(down_three_sigma_trend).mean()
   
# P( }/ d6 ], Z+ F% x1 Y" m/ K    return sims,midline_trend,up_one_sigma_trend,up_two_sigma_trend,up_three_sigma_trend,down_one_sigma_trend,down_two_sigma_trend,down_three_sigma_trendsims,midline_trend,up_one_sigma_trend,up_two_sigma_trend,up_three_sigma_trend,down_one_sigma_trend,down_two_sigma_trend,down_three_sigma_trend = simulation(model.decode(features)[1][-1],_states,10000,100,trans_mat,means_mat,covar_mat)
test_data = test_data.reset_index(drop=True)
frac_change_test = _extract_features(test_data)[1][:100]
frac_high_test = _extract_features(test_data)[2][:100]
frac_low_test = _extract_features(test_data)[3][:100]
3 j+ L% G( v1 ~+ v
( n$ o' Z& E5 Cdef get_values(arr):
    result = []
    for i in range(len(arr)):
        if i == 0:
* q0 ~: E2 y1 k  @- a  I: B            result.append(1*(1+arr[0]))
: n2 b( o% X) o5 S        else:
3 v& S' w; l  M            result.append(result[i-1]*(1+arr))
    return result

plt.figure(figsize=(15,10))
0 }. o! n+ [8 k- M5 [plt.plot_date(test_data[&#39;开盘时间_&#39;][:100].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;)),frac_change_test,&#39;-&#39;,label=&#39;close&#39;,lw=1)
plt.plot_date(test_data[&#39;开盘时间_&#39;][:100].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;)),frac_high_test,&#39;-&#39;,label=&#39;high&#39;,lw=1)
" `- e. l1 C& q+ f3 Q' r# r7 n# jplt.plot_date(test_data[&#39;开盘时间_&#39;][:100].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;)),frac_low_test,&#39;-&#39;,label=&#39;low&#39;,lw=1)
plt.plot_date(test_data[&#39;开盘时间_&#39;][:100].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;)),midline_trend,&#39;-&#39;,label=&#39;midline&#39;,lw=0.5)
plt.plot_date(test_data[&#39;开盘时间_&#39;][:100].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;)),up_one_sigma_trend,&#39;-&#39;,label=&#39;up-1-sigma&#39;,lw=0.5)
plt.plot_date(test_data[&#39;开盘时间_&#39;][:100].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;)),up_two_sigma_trend,&#39;-&#39;,label=&#39;up-2-sigma&#39;,lw=0.5)
plt.plot_date(test_data[&#39;开盘时间_&#39;][:100].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;)),up_three_sigma_trend,&#39;-&#39;,label=&#39;up-3-sigma&#39;,lw=0.5)
plt.plot_date(test_data[&#39;开盘时间_&#39;][:100].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;)),down_one_sigma_trend,&#39;-&#39;,label=&#39;down-1-sigma&#39;,lw=0.5)
plt.plot_date(test_data[&#39;开盘时间_&#39;][:100].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;)),down_two_sigma_trend,&#39;-&#39;,label=&#39;down-2-sigma&#39;,lw=0.5)
4 H! W, B" v% \! x7 h# l7 Jplt.plot_date(test_data[&#39;开盘时间_&#39;][:100].apply(lambda x:datetime.datetime.strptime(x, &#34;%Y-%m-%d %H:%M:%S&#34;)),down_three_sigma_trend,&#39;-&#39;,label=&#39;down-3-sigma&#39;,lw=0.5)
plt.legend()结果如下：

3 w* I1 S3 }$ U9 m- A3 I
) @2 Z" c6 g; @% O* \9 D; E 可以看出来大部分5min波动都在两个sigma之内，不过关于实际应用还没有想到有什么可以比较好利用这个点的方式，可以作为未来探索方向。  
) ~/ z2 g  B# ^" p9 z<hr/>6. 小结和Future Works

: ]; Z) l. o# N# e3 s6 ]+ F
虽然通过这种方式可以从另一个角度解读价格走势的变化，但是实际应用上还有所欠缺，目前做的这些研究工作还是比较碎片化，下次再做HMM的探索时希望可以沉淀一些策略层面上的东西。
一些Future Works：

7 ]) U1 M& t$ Z" I% C
• 使用价格的对数变化率替代绝对变化率；
  $ _! r& ?% G( T& @0 u7 ~! C
• 显性状态抽样优化，可以使用拟合后的分布进行抽样；
  , T% t$ c' U, ?6 _- \# S
• 使用更丰富的隐性状态来解释价格走势状态；
  
• 探索隐性状态转移在策略层面上的应用。
  " y) M7 a- `: W0 o

此外，对于HMM的应用层面如果有任何想法欢迎留言交流～  
谢谢大嘎！
% R2 g# D! {# r' j5 }附：前几次的文章
加密货币衍生品实战日记（八）—— 币本位溢价与资费套利策略详解
: N5 X/ h6 p1 ?6 M& z加密货币衍生品实战日记（七）——蒙特卡洛在价格走势范围预测上的应用(译)
" U# x: |" u& B$ y6 W8 X2 F. M& I+ _加密货币衍生品实战日记（六）——隐马尔科夫链对于状态的预测
<hr/>
* j( x, N2 `( s* W" g