python7c 高等数学 之 BL模型，barra收益回归模型

1. BL

def blacklitterman(returns,tau,P,Q):

mu=returns.mean()

sigma=returns.cov() #协方差

pil=np.expand_dims(mu,axis=0).T #转换维度，axis=0代表扩展列； axis=-1代表扩展最后一个参数

ts=tau*sigma

ts_1=linalg.inv(ts) #np.linalg.inv()：矩阵求逆

Omega = np.dot(np.dot(P,ts), P.T)* np.eye(Q.shape[0])#如果a和b都是一维数组，np.dot 向量的内积；如果a和b都是二维数组，那么它返回的是矩阵乘法。# np.eye 对角线为1，其他为0

Omega_1 = linalg.inv(Omega)

er = np.dot(linalg.inv(ts_1 + np.dot(np.dot(P.T,Omega_1),P)),(np.dot(ts_1 ,pil)+np.dot(np.dot(P.T,Omega_1),Q)))

posterirorSigma = linalg.inv(ts_1 + np.dot(np.dot(P.T,Omega_1),P))

return [er, posterirorSigma]

##定义最小化方差的函数,即求解二次规划

def blminVar(blres, goalRet):

covs = np.array(blres[1],dtype=float)

means = np.array(blres[0],dtype=float)

L1 = np.append(np.append(covs.swapaxes(0,1),[means.flatten()],axis=0),[np.ones(len(means))],axis=0).swapaxes(0,1)

## 使用transpose（）函数转换各个元素的维度，其括号内如果不加入参数，其含义为：将数组的各个元素，从原默认维度顺序，转换至括号内维度顺序。

##swapaxes类似，是作用在维度上的

# #means.flatten()， 返回一维数组

L2 = list(np.ones(len(means)))

L2.extend([0,0])

##数据拼接

L3 = list(means)

L3.extend([0,0])

L4 = np.array([L2,L3],dtype=float)

L = np.append(L1,L4,axis=0)

results = linalg.solve(L,np.append(np.zeros(len(means)),[1,goalRet]))

## b。a是一个N*N的二维数组，而b是一个长度为N的一维数组，solve函数找到一个长度为N的一维数组x，使得a和x的矩阵乘积正好等于b，数组x就是多元一次方程组的解。

# return pd.DataFrame(results[:-2],columns = ['p_weight'])

#blresult = blminVar(res,0.70/252)

#print(blresult)

对权重进行赋值

def statistics(weights):

weights = np.array(weights)

# returns.mean() = np.array(blres[1],dtype=float)

# returns.cov()= np.array(blres[0],dtype=float)

port_returns = np.sum(returns.mean()*weights)*252

port_variance = np.sqrt(np.dot(weights.T, np.dot(returns.cov()*252,weights)))

sharpe=port_returns/port_variance

return np.array([port_returns, port_variance, sharpe])

#最小化夏普指数的负值,port_returns/port_variance为夏普

def min_sharpe(weights):

return -statistics(weights)[2]

def con(args):

#约束是所有参数(权重)的总和为1。这可以用minimize函数的约定表达如下,type可选'eq'和'ineq'，分别是等式=0约束和不等式约束>=0，

x1min,x1max,x2min,x2max,x3min,x3max=args

cons = ({'type':'eq', 'fun':lambda x: np.sum(x)-1},\

{'type': 'ineq', 'fun': lambda x: x[0] - x1min},\

{'type': 'ineq', 'fun': lambda x: -x[0] + x1max},\

{'type': 'ineq', 'fun': lambda x: x[1] - x2min},\

{'type': 'ineq', 'fun': lambda x: -x[1] + x2max},\

{'type': 'ineq', 'fun': lambda x: x[2] - x3min},\

{'type': 'ineq', 'fun': lambda x: -x[2] + x3max})

return cons

args1 = (0.1,0.9,0.1, 0.9,0.1,0.9) #x1min, x1max, x2min, x2max

cons = con(args1)

#设置初始猜测值 ,

x0 = np.asarray((0.5,0.5,0.5))

res = sco.minimize(min_sharpe(weights), x0, method='SLSQP',constraints=cons)

print(res.fun)

print(res.success)

print(res.x)

2. Barra 回归

（1） 需要按昨日和今日计算的部分

（2） 去极值、标准化

（3） WLS 加权回归

barra_returns = pd.DataFrame()

for i in last_day.index: # 这里的last_day 是前一日， this_day 是当日

try: # 谨防barra日期没有及时更新

# 找到回归的x和y，收益率为模拟法下股票收益

# x 为前一日实际barra的因子值

c1, barra_raw1 = _get_temp_combine(c, last_day[i])

factor_c1 = factor_c1.set_index('code')

# factor_market 是 barra在前一日的所有股票因子值

factor_market = barra_raw1.set_index('code')

# pct_thisday 是当日的股票收益矩阵

pct_thisday = stock_pct.loc[this_day[i], :]

pct_thisday.index = [i.split(".")[0] for i in pct_thisday.index]

# mv是股票的市值矩阵

mv_thisday = stock_mv.loc[this_day[i], :]

mv_thisday.index = [i.split(".")[0] for i in mv_thisday.index]

# 去极值

returns_pct_market2 = pct_thisday * 0.01

def extreme_percentile(rawdata, min=0.025, max=0.975):

p = rawdata.quantile([min, max]) # 得到上下限的值

return rawdata.clip(p.loc[min, :], p.loc[max, :], axis=1) # 超出上下限的值，赋值为上下限

factor_market2a = extreme_percentile(factor_market, min=0.1, max=0.9)

# 标准化

factor_market2b = (factor_market2a.T - (factor_market2a.dropna()).mean(axis=1)) / (factor_market2a.dropna()).std(axis=1)

factor_market2 = factor_market2b.T.dropna()

# 直接使用weight进行加权， barra中已验证使用根号市值，进行加权

weight = (np.sqrt(mv_thisday) / (np.sqrt(mv_thisday)).sum())

weight1 = weight[factor_market2.index].replace(np.nan, 0).replace(-np.inf, 0).replace(np.inf, 0)

dataall_market2 =pd.merge(returns_pct_market2 , factor_market2, left_index=True, right_index=True, how='inner') # 将回归的股票数量对齐

dataall_market2 = dataall_market2.replace(np.nan, 0).replace(-np.inf, 0).replace(np.inf, 0)

factor0_market2 = dataall_market2.iloc[:, 1:] # 全市场barra因子，排除国家

returns0_market2 = dataall_market2.iloc[:, 0] # 全市场收益率

x = factor0_market2

y = returns0_market2

model = sm.WLS(y, x, weights=weight1).fit()

z = model.params[:]

barra_returns.loc[:, this_day[i]] = z # 需要按昨日和今日计算的部分