关于量化。－－镭矿

一、参数优化的 Bayesian Optimization 方法 未方便查看，可看原文链接 http://raquant.com/qa/index.php?qa=248&qa_1= 很多时候我们比较头痛参数该怎么选取的问题，比如： KDJ 的回看时间设多少， PE 的上限设多少，版块的选股比例设多少。即便不是处女座，在我们写交易策略的时候也都希望能把参数调整到一个自己满意的数值，当交易系统变的复杂，手动优化这些参数就显得非常的费力，而单用历史数据去一个模块一个模块的优化，很有可能会过度拟合。

那该如何是好呢？

这里我来简单分享下我自己的参数优化方法： Bayesian Optimization

先说明一下， Bayesian Optimization 和传统的优化方法完全不是一回事，一般我们优化就是梯度下降， variational inference 等解析方法或者 MCMC ， particle filter 等随机取样的方法。但这些方法的前提是这个函数你是知道的。在很多实际情况下，原函数是没有数学表达式的，比如交易系统的回测效果，他也是一个函数，但是没有表达式。 这个时候就需要 BO 了。

Bayesian Optimization ：

    Beyesian Optimization(BO)是在原函数未知的情况下去最优化该函数的一种方法。比如，我们有个交易策略，里面有三个参数需要优化来使得该策略回测的收益最大。那这里的函数就是我们的交易策略，输入是这三个参数，输出是回测的收益率。我们不可能用梯度方法找参数的最优解，因为交易策略木有梯度可求，臣妾算不出来啊。这时候就需要 BO 闪亮登场了。

   BO 是先用一个 data driven 的模型去近似原函数，然后在近似的函数上找最值。在回测的例子中， BO 会用一个模型去构建交易策略的参数和回测效果之间的函数关系，我们用不同的参数进行回测，会有不同的回测表现，再根据这些表现去调整参数。是不是和我们手动去猜参数的过程一模一样！！！

    对， BO 其实就是在做我们平常一直做的事情，不停的回测不停的猜。但是我们总希望通过尽可能少的回测次数来得到最优的参数。这时候盲目去猜就比 BO 的效果差距十万八千里了。这里我就简单介绍下神奇而贴心的 BO 。

import numpy as np import pandas as pd import matplotlib.pyplot as plt import scipy as sp from sklearn import gaussian_process from raquant import * import math 假设我们有如下一个未知的函数，我们不知道其函数形式，只能是给出自变量 xx ，然后函数能返回一个应变量 yy ，其他信息一概不知。

为了实验起见，我们假装不知道如下一个函数：

def func(x): mu1 = 0.3 var1 = 0.01 mu2 = 0.6 var2 = 0.03 y = (1/np.sqrt(2math.pivar1)) * math.exp(-(x-mu1)2./var1) + (1/np.sqrt(2math.pivar2)) * math.exp(-(x-mu2)**2./var2) return y 假设我们只知道如下的 3 个样本点了：

x = np.linspace(0,1,200) y=[] for i in x: y.append(func(i)) y=np.array(y) plt.plot(x,y,'-r')

xtr=np.array([0.1,0.5,0.8]) ytr=[] for i in xtr: ytr.append(func(i)) ytr=np.array(ytr) plt.plot(xtr,ytr,'ko',markersize = 10)

当前只有三个观测值 (x=0.1 ， y=0.0736),(x=0.5,y=1.7235),(x=0.8,y=0.6071)(x=0.1 ， y=0.0736),(x=0.5,y=1.7235),(x=0.8,y=0.6071)

我们的目标是尝试不同的 xx 来找到 yy 的最小值。再强调下，如果我们知道函数的形式，可以用梯度下降方法直接解决，但问题是我们连函数长啥样都没见过，只能瞎猫碰死老鼠了。那怎么去猜，如何减少猜的次数就是个问题。

好，既然函数本身不知道，我们用其他的函数模拟他， BO 常用的拟合工具就是万金油的 Gaussian process(GP)。为啥用 GP 呢，原因很简单，在我们对建模缺乏经验时都会无脑假设为高斯类分布。。。现在做 BO 的也有拿 deep network 来做的，不过纯属因为 deep learning 比较火，好发文章，其实实际效果不如 GP 来的靠谱。

直接上公式吧，给定观测值 yD,对某个给定的 xx ，其观测值 yxyx 的后验分布 N(yx|D,Kxx|D)如下：

yx|D=μx+KxDK−1DD(y−μD)

Kxx|D=Kxx−KxDK−1DDKDx

可能你已经注意到了，我们其实不需要把函数完全拟合出来，我们只需要把可能出现最大值的区域拟合好就行了，为了达到这个效果， BO 需要处理一个经典的问题 exploration vs. exploitation 。意思是说，我现在应该去探索未知的领域，还是最优化我当前的结果。我今天该去吃没吃过的菜，还是吃我最爱吃的菜。。。

BO 使用 acquisition function 来处理这个问题。比较常见的 acquisition function 有 expected improvement 和 UCB 。我个人偏爱 UCB ，因为相关的证明比较好推导。 UCB 的形式是：

UCB(x)=μx|D+βKxx|

BO 每次取 UCB(x)UCB(x)最大的那个 xx 来进行实验，这样可以保证函数的拟合度，又能尽快找到最优解。至于 UCB 为啥能达到这个效果，就需要看相关文献了，他是通过ββ的不同取值办法，得到一个 probability bound ，然后这个 bound 可以保证 BO 整个过程可以收敛。我们在此只是应用 BO ，所以不必太纠结这个理论推导。

我们就拿刚才的三个观测值来作为 GP 的输入，看看计算出来的函数长啥样。其中，红线是真实的函数，蓝色线是我们预测的函数值，绿色的线是 UCB 的值。

D = np.atleast_2d(xtr).T y_D = [] for i in xtr: y_D.append(func(i)) y_D=np.array(y_D) gp = gaussian_process.GaussianProcess(theta0=0.05, thetaL=0.01, thetaU=1) gp.fit(D, y_D)

xt = np.atleast_2d(x).T y_pred, sigma2_pred = gp.predict(xt, eval_MSE=True)

这个红线和蓝线看着差别略大啊，完全不是一个妈生的，有木有！！！不过当我们有更多的观测值的时候，效果会慢慢出来的。但是需要强调一点，那个绿线，也就是 UCB(x)，在观测值附近会比较接近蓝线，在没有观测值的地方会偏离蓝线。这个偏离的程度也就是我们对原函数的把握程度，偏离的越小，说明我们拟合的越好。而在最值的附近， UCB(x)会有比较大的值，上图也是这样的情形，所以我们只要每次选绿线最大的点做实验就会慢慢逼近最优解了。

xtr2 = np.random.sample(40) D = np.atleast_2d(xtr2).T y_D = [] for i in xtr2: y_D.append(func(i)) y_D=np.array(y_D)

gp = gaussian_process.GaussianProcess(theta0=0.05, thetaL=0.01, thetaU=1) gp.fit(D, y_D)

xt = np.atleast_2d(x).T y_pred, sigma2_pred = gp.predict(xt, eval_MSE=True)

plt.plot(x,y,'-r') plt.plot(xt,y_pred,'-b') plt.plot(xt,y_pred+2*np.sqrt(sigma2_pred),'-g')

当我们取了 40 个观测点，效果就出来了，简直就是双胞胎，完全分不出来了。当观测值原来越多时，虽然我们不知道原函数，但 GP 可以基本把原函数给描述出来，我们只需要找这个 GP 的最大值就好了。那个绿线 UCB(x)现在基本和蓝线完全贴合，也就是说我们对原函数的拟合非常好了，同时绿线最大值的点就是原函数最优解的点。

在回测中应用 BO 时，可以这么做：先自己随便设几个参数值跑几遍回测，把观测结果输入到 BO 中，他会给出下一次回测需要设置的参数值，你再手动去跑一次回测，再把回测结果返回给 BO,如此往返几百次，就能得到最优解。一般一个参数需要 50 次的实验，两个参数需要 200 次实验。。。我们交易策略里的参数通常不下 10 个，如果完全手动去这么做，臣妾做不到啊。

二、 hurst 指数 为方便查看，可看原文链接 http://raquant.com/qa/index.php?qa=192&qa_1=hurst 本文将对 hurst 指数的历史及数学背景做尽可能详细而直白的介绍。

简介：

1.布朗运动（ Brownian motion ）：被分子撞击的悬浮微粒做无规则运动的现象叫做布朗运动。例如常见的解释是：显微镜下水中花粉微粒在水分子的撞击下观测到的运动轨迹。布朗运动是大量分子做无规则运动对悬浮的固体微粒各个方向撞击作用的不均衡性造成的，所以布朗运动是大量液体分子集体行为的结果，是一种随机涨落现象。由布朗发现，后经爱因斯坦、维纳等人发展成为分子运动论和统计力学的基础。

2.布朗运动数学描述：

{B(t)}布朗运动，满足下面性质：

  1.B(t)是独立增量过程：对于任意的 t>s, B(t)-B(s)独立于之前的过程 B(u):0<=u<=s.

  2. B(t)有正态的增量： B(t)-B(s)满足均值为 0 方差为 t-s 的正态分布。

  3.B(t)有连续的路径： B(t), t>=0 是关于 t 的连续函数。

下图是二维布朗运动的一个模拟图，像不像像花粉微粒在水中的运动咯~wink

我们再看看一维的标准布朗运动：均值为 0 ，方差为 sqrt （ T ） ,（你能找到它和股票价格走势的区别我给你五毛钱~laugh ）

3.爱因斯坦的发现

爱因斯坦于 1908 年发表了关于布朗运动的论文，使得布朗运动成为随机游走的基本模型。爱因斯坦发现，分子布朗运动的半径可以用涵盖时间的平方根来测度：

                         R=T^(0.5)；                 其中 R 分子随机运动的半径（如上图）， T 是时间（ Time ）     （ 1 ）

Hurst 指数就是在这基础上发现新的统计量来对时间序列的随机性进行检验。

4.Hurst 指数背景

由英国水文学家 H.E.HURST （赫斯特）提出，用于研究洪水形成周期：洪水过程是时间系列曲线，具有正的长时间相关效应。即干旱愈久，就可能出现持续的干旱；大洪水年过后仍然会有较大洪水。（有木有发现和股市很像 cool ？！）

而 Hurst 发现上述这种特性可以用他本人提出的 Hurst index 来衡量。

5.Hurst 指数数学依据

从上述爱因斯坦的结论（ 1 ）式中，我们知道标准布朗运动的半径 R=T^(0.5)，那么对于一般的时间序列呢？！受爱因斯坦的结果影响，我们相信分子运动半径 R 与时间 T 依旧存在类似的关系： R~T^(H)，其中 H 就是所谓的 Hurst 指数，这里不再用等号是认为此时半径应该正比于 T^(H)。

经研究得出 Hurst 指数的三种形式：

1 ．如果 H=0.5 ：标准布朗运动，表明时间序列可以用随机游走来描述，表现出马尔科夫性；

2 ．如果 0.5<H<1 ，表明时间序列存在长期记忆性，表明一定时间范围内走势记录会持续相当长的时间，从而形成一个大的循环；

3 ．如果 0≤H<0.5 ，表明粉红噪声(反持续性)即均值回复过程，该序列理论上无数次返回它的历史出发点（说的就是此时股票价格就是处在震荡时刻）。

也就是说，只要 H≠0.5 ，就可以用有偏的布朗运动(分形布朗运动)来描述该时间序列数据。

问题来了，那怎么计算 Hurst 指数呢？（看完下面的计算过程，我想你会后悔这么着急着想看这了~angry ，哈，逗你的，计算根本不难，根本不难，根本不难~如果你耐心看的话！）

6.Hurst 指数计算过程（以股价时间序列为例）

7.这样就完成 Hurst index 的计算了咯~哎，别走啊，再看一遍吧，真的不难啊~broken heart

下面我以上证指数为例，计算了 2013 年以来该指数 N （=120 ）日的 Hurst 指数：

我们可以看到该上证指数大概在 2/3 的时间里的 Hurst 指数是大于 0.5 的，也就是说，上证指数的走势在很大范围里具有长周期记忆过程，而在约 1/3 的时刻处于震荡状态~（这可不是我说的~不信你看图 wink ）

8.下面是该图的源代码：

import numpy as np import pandas as pd from scipy import stats from statsmodels import regression import math

def init(context): #stocks=find_by_group("sz50") stocks="sha-6000001" universe.extend(stocks) context.bench_mark="sza-000001" context.slippage=0.003 context.N=121 #取样数 context.A=range(2,61)

def handle_data(context,data): stock=context.bench_mark n=context.N his=history(n,"close")[stock] his1=[] for i in range(0,len(his)-1): his1.append(np.log(float(his[i+1])/float(his[i]))) x=[] y=[] for a in context.A: M=(n-1)/a x.append(np.log(a)) R_S=0 for i in range(M): his2=his1[i*a:(i+1)*a] mean=np.mean(his2) st=np.sqrt(np.var(his2)) his3=[0]*a his3[0]=his2[0]-mean for j in range(1,a): his3[j]=his3[j-1]+his2[j]-mean R=max(his3)-min(his3) R_S=R_S+R/st y.append(np.log(R_S/M)) slope, intercept, r_value, p_value, slope_std_error = stats.linregress(x, y) H=slope #log.info("x===%s"%x) record("H",H) #log.info("y===%s"%y) #log.info("H==%s"%H) #log.info("%s"%data.dt)

有了上面这个计算过程，至于怎么形成交易信号就比较容易了咯~（要不你先试试？！ yes ）

这里我就用 sha-600000 作为参考，以 H>0.62 作为买入信号， H<0.48 作为卖出信号。可以看到建仓出仓次数很少，至于效果嘛~我没有过多尝试尝试的设置，自己动手做做咯~这里只是参考嘛，当然你也可以重新设置股票或者以股票池来操作，具体这些都很简单啦~

源代码及回测：

import numpy as np import pandas as pd from scipy import stats from statsmodels import regression import math

def init(context): stocks="sha-600000" universe.append(stocks) context.bench_mark="sha-600000" context.slippage=0.003 context.N=121 #取样数 context.A=range(2,61)

def handle_data(context,data): stock="sha-600000" n=context.N his=history(n,"close")[stock] his1=[] for i in range(0,len(his)-1): his1.append(np.log(float(his[i+1])/float(his[i]))) x=[] y=[] for a in context.A: M=(n-1)/a x.append(np.log(a)) R_S=0 for i in range(M): his2=his1[i*a:(i+1)*a] mean=np.mean(his2) st=np.sqrt(np.var(his2)) his3=[0]*a his3[0]=his2[0]-mean for j in range(1,a): his3[j]=his3[j-1]+his2[j]-mean R=max(his3)-min(his3) R_S=R_S+R/st y.append(np.log(R_S/M)) slope, intercept, r_value, p_value, slope_std_error = stats.linregress(x, y) H=slope #log.info("x===%s"%x) record("H",H) #log.info("y===%s"%y) #log.info("H==%s"%H) #log.info("%s"%data.dt) if H>0.62: if context.portfolio.positions[stock].amount==0: cash=context.portfolio.cash order_value(stock,cash) log.info("Buying~") if H<0.48: if context.portfolio.positions[stock].amount>0: order_target(stock,0) log.info("selling~")

期待各路豪杰疯狂回复