【人工智能与量化交易】浅谈贝叶斯统计 | RiceQuant米筐量化社区 交易策略论坛

前言篇

理论篇

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  #记录出现正面的次数
count_0 = 0  #记录出现反面的次数
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)
         #画出概率密度函数
        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   #出现正面
    else:
        count_0 += 1   #出现反面
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()

实战篇

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

#现在假设我们想估计近期以下4个股票的期望收益率
stock_dict = {'贵州茅台':'600519.SH',
              '农业银行':'601288.SH',
              '中国石油':'601857.SH',
              '工商银行':'601398.SH'}
daily_log_return = {}
#从tushare上获取数据
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')
    #每日的对数收益率
    daily_log_return[stock] =  np.log(df['close']/df['pre_close'])  
daily_log_return = pd.DataFrame(daily_log_return).sort_index()
#设置各个股票收益率期望值所服从的先验分布，这里仅作为例子示范
export_prior_params = stock_dict = {'贵州茅台':(0.004,0.02),
                                    '农业银行':(0.005,0.02),
                                    '中国石油':(0.006,0.01),
                                    '工商银行':(0.003,0.015)}
#调整好dataframe列的顺序
daily_log_return = daily_log_return[[x for x in 
                                     export_prior_params.keys()]]
#先验分布的均值和标准差
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
#从正态分布中生成均值向量
mu = pm.Normal('returns',prior_mu,1,size=4)
#从威沙特分布中生成协方差矩阵的逆矩阵
inv_cov_matrix = pm.Wishart('inv_cov_matrix',
                            len(daily_log_return),
                            np.diag(prior_std**2)) 
#由均值向量和协方差矩阵的逆生成观测值
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组参数并从10000组以后开始采样
#对样本求均值,得到参数的估计
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):
    '''
    协方差矩阵转化为相关系数矩阵
    '''
    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')