import numpy as np
import pandas as pd
from arch import arch_model # GARCH(1,1)
from matplotlib import pyplot as plt
from datetime import timedelta
from CAL.PyCAL import *
from lib.Hurst import * inter = 320 #滑动时间窗口
#设置时间
today = Date.todaysDate()
beginDate = '20100101'
endDate = today.toDateTime().strftime('%Y%m%d') #设置指数类型
indexLabel = '000001' # SSE index
#indexLabel = '399006' # CYB index#读取指数
indexPrice = DataAPI.MktIdxdGet(ticker=indexLabel,beginDate=beginDate,endDate=endDate,field=["tradeDate","closeIndex"],pandas="1")
price = np.array(indexPrice.loc[:,'closeIndex'])#计算对数收益
back_price = np.append(price[0],price.copy())
back_price = back_price[:-1]
return_price = np.log(price) - np.log(back_price)#计算波动率 from GARCH(1,1)
am = arch_model(return_price)
res = am.fit()
sqt_h = res.conditional_volatility#去除波动性
f = return_price/sqt_h#计算hurst指数,函数来自自定义library
hurst = Hurst(f,T=inter,step=1,q=2,Smin=10,Smax=50,Sintr=1)indexPrice['Hurst'] = pd.DataFrame(np.array([0] * len(indexPrice)))
indexPrice.loc[inter-1:,'Hurst'] = hurst
indexPrice.index = indexPrice['tradeDate']

plt.figure(figsize=(10,6))
plt.subplot(3,1,1)
plt.plot(f)
plt.subplot(3,1,2)
plt.plot(return_price)
plt.subplot(3,1,3)
plt.plot(sqt_h)

plt.figure(1)
indexPrice['closeIndex'].tail(len(indexPrice)-inter).plot(figsize=(10,4),color='red',title='SSE Index',linewidth=1)
plt.figure(2)
indexPrice['Hurst'].tail(len(indexPrice)-inter).plot(figsize=(10,4),color='green',title='Hurst Index',linewidth=1,marker='.')

• 看出了啥没？简单点说，就是hurst越大，越有可能延续之前的趋势（即动量），若hurst越小，则越有可能违反之前的趋势（即反转）。LADISLAV KRISTOUFEK这位教授的想法是通过极大极小值来判断，当然它分析的是美股啦。
• 再看看上面的图，是对上证指数的分析，取的是日线的数据（其实我喜欢用分钟线，因为A股波动辣么牛逼，日线颗粒度哪里够啊。。），可以得（meng）出这些结论：
  

  1. 13年中旬hurst出现最小值，说明熊市的跌势要反转了，马上要进入牛市了？！
  2. 15年中旬hurst出现最小值，说明牛市的涨势要反转了，马上要进入熊市了？！

• 算卦完毕。

#coding=utf-8import numpy as np
from sklearn.linear_model import LinearRegressiondef avgFluctuation(Xseries,q,S):lr = LinearRegression(fit_intercept = True)T = Xseries.shape[0]Ts = int(T/S)qorder = qaFlu = 0Xtime = np.array(range(1,S+1))Xtime.shape = (S,1)for v in range(1,Ts+1):Xstarts = (v-1)*SXends = v*SXtemp = Xseries[Xstarts:Xends]Xtemp.shape = (S,1)lr.fit(Xtime,Xtemp)fitX = map(lr.predict,Xtime)eX = fitX - XtempFlu = np.sum(eX*eX)/SaFlu += (Flu)**(qorder/2)/(2*Ts)for v in range(Ts+1,2*Ts+1):Xstarts = T-(v-Ts)*SXends = T-(v-Ts)*S+SXtemp = Xseries[Xstarts:Xends]Xtemp.shape = (S,1)lr.fit(Xtime,Xtemp)fitX = map(lr.predict,Xtime)eX = fitX - XtempFlu = np.sum(eX*eX)/SaFlu += (Flu)**(qorder/2)/(2*Ts)aFlu = aFlu**(1./qorder)return aFludef MFDFA(Xseries,q,Smin,Smax,Sintr=1):T = Xseries.shape[0]N = (Smax - Smin + 1)/SintraFlus = np.zeros(N)i = 0for S in range(Smin,Smax+1,Sintr):aFlus[i] = avgFluctuation(Xseries,q,S)i += 1logaFlus = np.log(aFlus)logS = np.log(range(Smin,Smax+1,Sintr))lr = LinearRegression(fit_intercept = True)logaFlus.shape = (N,1)logS.shape = (N,1)lr.fit(logS,logaFlus)h = lr.coef_return hdef Hurst(X,T=300,step=1,q=2,Smin=10,Smax=50,Sintr=1):X = np.array(X)nX = X.shape[0]hurst = np.zeros(nX-T+1)for i in range(0,nX-T+1,step):XX = X[i:i+T]Xseries = np.zeros(T)for j in range(T):Xseries[j] = np.sum(XX[0:j+1])hurst[i] = MFDFA(Xseries,q,Smin,Smax,Sintr)return hurst 

def avgFluctuation(Xseries, q, S):lr = LinearRegression(fit_intercept = True)T = Xseries.shape[0]Ts = int(T/S)qorder = qaFlu = 0Xtime = np.array(range(1, S+1))Xtime.shape = (S, 1)for v in range(1, Ts+1):Xstarts = (v-1)*SXends = v*SXtemp = Xseries[Xstarts:Xends]Xtemp.shape = (S, 1)lr.fit(Xtime, Xtemp)fitX = list(map(lr.predict, [Xtime]))eX = fitX - XtempFlu = np.sum(eX*eX) / SaFlu += (Flu)**(qorder/2) / (2*Ts)for v in range(Ts+1,2*Ts+1):Xstarts = T-(v-Ts)*SXends = T-(v-Ts)*S+SXtemp = Xseries[Xstarts:Xends]Xtemp.shape = (S,1)lr.fit(Xtime,Xtemp)fitX = list(map(lr.predict, [Xtime]))eX = fitX - XtempFlu = np.sum(eX*eX) / SaFlu += (Flu)**(qorder/2) / (2*Ts)aFlu = aFlu**(1./qorder)return aFludef MFDFA(Xseries, q, Smin, Smax, Sintr=1):T = Xseries.shape[0]N = int((Smax - Smin + 1) / Sintr)aFlus = np.zeros(N)i = 0for S in range(Smin, Smax+1, Sintr):aFlus[i] = avgFluctuation(Xseries, q, S)i += 1logaFlus = np.log(aFlus)logS = np.log(range(Smin, Smax+1, Sintr))lr = LinearRegression(fit_intercept = True)logaFlus.shape = (N,1)logS.shape = (N,1)lr.fit(logS, logaFlus)h = lr.coef_return hdef Hurst(X, T=300, step=1, q=2, Smin=10, Smax=50, Sintr=1):X = np.array(X)nX = X.shape[0]hurst = np.zeros(nX-T+1)for i in range(0, nX-T+1, step):XX = X[i:i+T]Xseries = np.zeros(T)for j in range(T):Xseries[j] = np.sum(XX[0:j+1])hurst[i] = MFDFA(Xseries, q, Smin, Smax, Sintr)return hurst