【人工智能與量化交易】淺談貝葉斯統(tǒng)計 | RiceQuant米筐量化社區(qū) 交易策略論壇

烏商匯富 2019-05-22

展開全文

前言篇

理論篇

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
import scipy.stats

p = 0.5  #上帝的聲音
count_1 = 0  #記錄出現(xiàn)正面的次數(shù)
count_0 = 0  #記錄出現(xiàn)反面的次數(shù)
plot_n = [0,1,2,3,4,5,8,15,50,500]
x = np.arange(0,1.01,0.01)
plt.figure(figsize=(10,10))
np.random.seed(0)
for i in range(500+1):
    if i in plot_n:
        beta = scipy.stats.beta(a=count_1+1,b=count_0+1)   #Beta分布
        plt.subplot(5,2,plot_n.index(i)+1)
        y = beta.pdf(x)
         #畫出概率密度函數(shù)
        plt.plot(x,y,color='blue',
                 label='觀測{}次,得到{}次正面'.format(i,count_1))
        plt.fill_between(x,y,color='blue')
        plt.axvline(x=p,ymin=0,ymax=y.max()*1.5,color='red',ls='--')
        plt.ylim(0,y.max()*1.5)
        plt.legend(loc='upper right',fontsize='large')
        
    if np.random.rand()<p:
        count_1 += 1   #出現(xiàn)正面
    else:
        count_0 += 1   #出現(xiàn)反面
plt.show()

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
matrix = np.matrix([[0.9,0.075,0.025],
                    [0.15,0.8,0.05],
                    [0.25,0.25,0.5]], dtype=float)

vector1 = np.matrix([[1,0,0]], dtype=float)
vector2 = np.matrix([[0,1,0]], dtype=float)
vector3 = np.matrix([[0,0,1]], dtype=float)
proba1,proba2,proba3 = [],[],[]
for i in range(20):
    proba1.append(vector1)
    proba2.append(vector2)
    proba3.append(vector3)
    vector1 = vector1*matrix
    vector2 = vector2*matrix
    vector3 = vector3*matrix
proba = np.array([proba1,proba2,proba3]).reshape((3,20,3))
colors = ['red','green','blue']
state = ['bull','bear','stagnant']
t = np.arange(20,dtype=int)
plt.figure(figsize=(12,3))
for i in range(3):
    plt.subplot(1,3,int(i+1))
    for k in range(3):
        plt.plot(t,proba[i][:,k],color=colors[k],label=state[k])
    plt.ylabel('Probaility')
    plt.title('Start with '+state[i])
    plt.legend(loc='upper right')
plt.show()

實戰(zhàn)篇

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

#現(xiàn)在假設(shè)我們想估計近期以下4個股票的期望收益率
stock_dict = {'貴州茅臺':'600519.SH',
              '農(nóng)業(yè)銀行':'601288.SH',
              '中國石油':'601857.SH',
              '工商銀行':'601398.SH'}
daily_log_return = {}
#從tushare上獲取數(shù)據(jù)
import tushare as ts
pro = ts.pro_api('9c5d921fa5dd85b5a8ecf8c7b63092b5da313fcef98ff9185580c78e')
for stock in stock_dict.keys():
    df = pro.daily(ts_code=stock_dict[stock],
                   start_date='20180901', end_date='20180922')
    df= df.set_index('trade_date')
    #每日的對數(shù)收益率
    daily_log_return[stock] =  np.log(df['close']/df['pre_close'])  
daily_log_return = pd.DataFrame(daily_log_return).sort_index()
#設(shè)置各個股票收益率期望值所服從的先驗分布，這里僅作為例子示范
export_prior_params = stock_dict = {'貴州茅臺':(0.004,0.02),
                                    '農(nóng)業(yè)銀行':(0.005,0.02),
                                    '中國石油':(0.006,0.01),
                                    '工商銀行':(0.003,0.015)}
#調(diào)整好dataframe列的順序
daily_log_return = daily_log_return[[x for x in 
                                     export_prior_params.keys()]]
#先驗分布的均值和標(biāo)準(zhǔn)差
prior_mu = np.array([x[0] for x in export_prior_params.values()])
prior_std = np.array([x[1] for x in export_prior_params.values()])

#使用pymc庫
import pymc as pm
#從正態(tài)分布中生成均值向量
mu = pm.Normal('returns',prior_mu,1,size=4)
#從威沙特分布中生成協(xié)方差矩陣的逆矩陣
inv_cov_matrix = pm.Wishart('inv_cov_matrix',
                            len(daily_log_return),
                            np.diag(prior_std**2)) 
#由均值向量和協(xié)方差矩陣的逆生成觀測值
obs = pm.MvNormal('observed returns',
                  mu,inv_cov_matrix,
                  observed=True,value=daily_log_return.values)

model = pm.Model([obs,mu,inv_cov_matrix])
mcmc = pm.MCMC(model)
mcmc.sample(15000,10000) #生成15000組參數(shù)并從10000組以后開始采樣
#對樣本求均值,得到參數(shù)的估計
return_mean = mcmc.trace('returns')[:].mean(axis=0)
return_cov = np.linalg.inv(mcmc.trace('inv_cov_matrix')[:].mean(axis=0))
return_std = np.sqrt(return_cov.diagonal())



def cov2corr(A):
    '''
    協(xié)方差矩陣轉(zhuǎn)化為相關(guān)系數(shù)矩陣
    '''
    d = np.sqrt(A.diagonal())
    A = (A.T/d).T/d
    return A

import scipy.stats
plt.figure(figsize=(10,7))
for i in range(4):
    norm = scipy.stats.norm(return_mean[i],return_std[i])
    x = np.linspace(-0.05,0.05,100)
    y = norm.pdf(x)
    plt.subplot(2,2,int(i+1))
    plt.plot(x,y)
    plt.title(list(export_prior_params.keys())[i])
plt.show()

plt.imshow(cov2corr(return_cov),cmap=plt.cm.hot)
plt.xticks(np.arange(4),export_prior_params.keys())
plt.yticks(np.arange(4),export_prior_params.keys())
plt.colorbar(orientation='vertical')