#!/usr/bin/env python
# coding: utf-8

# ##内容索引
# 1. 相关性分析 --- cov函数、diagonal函数、trace函数、corrcoef函数
# 2. 多项式拟合 --- polyfit函数、polyval函数、roots函数、polyder函数
# 3. 计算净额成交量 --- sign函数、piecewise函数
# 4. 模拟交易过程 --- vectorize函数、round函数
# 5. 数据平滑 --- hanning函数

# In[2]:


get_ipython().run_line_magic('matplotlib', 'inline')
import numpy as np
from matplotlib.pyplot import plot
from matplotlib.pyplot import show


# ##1. 股票相关性分析
# 本例子中， 我们使用2个示例数据集提供收盘价数据，第一家公司是BHP Billiton(BHP)，其主要业务是石油、金属和钻石开采；第二家公司是Vale(VALE)，也是一家金属开采业公司，他们有部分业务是重合的。我们来分析一下他们的股票相关性。

# In[3]:


# 首先读入两只股票的收盘价，并计算收益率
bhp_cp = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
vale_cp = np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)
bhp_returns = np.diff(bhp_cp) / bhp_cp[:-1]
vale_returns = np.diff(vale_cp) / vale_cp[:-1]


# 协方差描述的是两个变量共同变化的趋势，其实就是归一化前的相关系数。
# 使用cov函数计算**股票收益率的协方差矩阵**

# In[4]:


covariance = np.cov(bhp_returns, vale_returns)
print 'Covariance:\n', covariance


# In[5]:


# 查看协方差矩阵对角线的元素
print 'Covariance diagonal:\n', covariance.diagonal()
# 计算矩阵的迹，即对角线之和
print 'Covariance trace:\n', covariance.trace()
# 计算相关系数，相关系数是协方差除以各自标准差的乘积
print 'Correlation coefficient:\n', covariance / (bhp_returns.std() * vale_returns.std())


# 用相关系数来度量两只股票的相关程度。相关系数的取值范围在-1到1之间，一组数据域自身的相关系数为1.使用corrcoef函数计算相关系数。

# In[6]:


# 使用corrcoef计算更加精确
print 'Correlation coefficient:\n', np.corrcoef(bhp_returns, vale_returns)


# 相关系数矩阵关于对角线对称，BHP与VALE的相关系数等于VALE和BHP的相关系数。看起来0.68的相关系数表示他们的相关程度似乎不是很强。

# ###判断两只股票的价格走势是否同步
# 如果它们的差值偏离了平均差值2倍于标准差的距离，则认为他们走势不同步

# In[7]:


difference = bhp_cp - vale_cp
avg = np.mean(difference)
dev = np.std(difference)
# 检查最后一次收盘价是否在同步状态
print "Out of sync : ", np.abs(difference[-1] - avg) > 2*dev


# 这说明，最后一次收盘价不再同步状态，我们暂时不能进行交易

# In[8]:


# 绘制收益率曲线
t = np.arange(len(bhp_returns))
plot(t, bhp_returns, lw=1)
plot(t, vale_returns, lw=2)
show()


# ##2. 多项式拟合
# NumPy中的plotfit函数可以用多项式去拟合一系列数据点，无论这些数据点是否来自连续函数都适用。

# In[9]:


# 用三次多项式去拟合两只股票收盘价的差价
t = np.arange(len(bhp_cp))
poly = np.polyfit(t, bhp_cp-vale_cp, 3)
print "Polynomial fit\n", poly


# In[10]:


# 用刚才得到的多项式对象，推断下一个值
print "Next value: ", np.polyval(poly, t[-1]+1)


# 理想情况下，BHP和VALE股票收盘价的差价越小越好。在极限情况下，差值可以在某个点为0。用roots函数找到拟合多项式函数在什么时候达到0。

# In[11]:


print "Roots: ", np.roots(poly)


# ###求极值

# In[12]:


# 极值位于导数为0的点
der = np.polyder(poly)
print "Dervative:\n", der
# 得到多项式导函数的系数


# In[13]:


# 求出导函数的根，即找出原多项式函数的极值点
print "Extremas: ", np.roots(der)


# In[14]:


# 通过argmax和argmin函数找到最大最小值点来检查结果
vals = np.polyval(poly, t)
print "Maximum index: ", np.argmax(vals)
print "Minimum index: ", np.argmin(vals)


# ###绘制拟合曲线

# In[15]:


plot(t, bhp_cp-vale_cp)
plot(t, vals)
show()


# ##3. 计算净额成交量
# 成交量表示价格波动的大小，净额成交量（On-Balance Volume)是由当日收盘价、前一天的收盘价以及当日成交量计算得出的。
# 
# 以前一日为基期计算当日的OBV的值（可认为基期的OBV的值为0）。若当日收盘价高于前一日收盘价，则本日OBV等于基期OBV加上当日成交量，否则减去当日成交量。

# **我们需要在成交量前面乘上一个有收盘价变化决定的正负号。**

# In[16]:


cp, volume = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,7), unpack=True)
change = np.diff(cp)
print "Change:", change


# 使用NumPy的sign函数返回每个元素的正负号。

# In[17]:


signs = np.sign(change)
print "Signs:\n", signs


# 使用Numpy的piecewise函数获取数组元素的正负。piecewise(分段)，可以根据给定取值，得到分段。

# In[18]:


pieces = np.piecewise(change, [change<0, change>0], [-1,1])
print "Pieces:\n", pieces


# In[19]:


# 检查两次输出是否一致
print "Arrays equal?", np.array_equal(signs, pieces)


# In[20]:


# OBV值的计算依赖于前一日的收盘价
print "On balance volume: \n", volume[1:]*signs


# ##4. 模拟交易过程

# 使用vectorize函数可以减少你的程序中使用循环的次数，NumPy中的vectorize函数相当于Python中的map函数。我们用它来计算单个交易日的利润。

# In[21]:


# 读入数据
# op is opening price,hp is the highest price
# lp is the lowest price, cp is closing price
op, hp, lp, cp = np.loadtxt('BHP.csv', delimiter=',', usecols=(3,4,5,6), unpack=True)


# 我们尝试以比开盘价稍低一点的价格买入股票。如果这个价格不在当日的股价范围内，则尝试买入失败，没有获利，也没有亏损，我们返回0。否则，我们将以当日收盘价卖出，所获得的利润即买入卖出的差价。

# In[22]:


def calc_profit(op, high, low, close):
    # 以开盘价买入,这里不考虑买入多少股
    buy = op
    
    if low < buy < high:
        return (close-buy) / buy
    else:
        return 0


# In[23]:


# 矢量化一个函数，这样可以避免使用循环
func = np.vectorize(calc_profit)
profits = func(op, hp, lp, cp)
print 'Profits:\n', profits


# In[24]:


# 我们选择非零利润的交易日并计算平均值
real_trades = profits[profits != 0]
print 'Number of trades:\n', len(real_trades), round(100.0 * len(real_trades)/len(cp), 2),"%"
print "Average profit/loss % :", round(np.mean(real_trades) * 100, 2)


# In[25]:


# 选择正盈利的交易日并计算平均利润
winning_trades = profits[profits > 0]
print "Number of winning trades", len(winning_trades), round(100.0*len(winning_trades)/len(cp),2),"%"
print "Average profit %", round(np.mean(winning_trades) *100, 2)


# In[26]:


# 选择负盈利的交易日并计算平均损失
losing_trades = profits[profits < 0]
print "Number of winning trades", len(losing_trades), round(100.0*len(losing_trades)/len(cp),2),"%"
print "Average profit %", round(np.mean(losing_trades) *100, 2)


# ##5. 数据平滑
# 噪声数据往往很难处理，我们需要对其进行平滑处理。

# hanning函数式一个加权余弦的窗函数，我们使用hanning函数平滑股票收益率的数组。
# 
# 离散卷积运算函数convolve的文档[convolve函数文档](http://docs.scipy.org/doc/numpy/reference/generated/numpy.convolve.html#numpy.convolve)

# In[27]:


# 调用hanning函数计算权重，生成一个长度为N的窗
# 这里N为8
N = 8
weights = np.hanning(N)
print "Weights:\n", weights


# In[28]:


# 首先读入两只股票的收盘价，并计算收益率
bhp_cp = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
vale_cp = np.loadtxt('VALE.csv', delimiter=',', usecols=(6,), unpack=True)
bhp_returns = np.diff(bhp_cp) / bhp_cp[:-1]
vale_returns = np.diff(vale_cp) / vale_cp[:-1]

# convolve函数，离散线性卷积运算
smooth_bhp = np.convolve(weights/weights.sum(), bhp_returns)[N-1 : -N+1]
smooth_vale = np.convolve(weights/weights.sum(), vale_returns)[N-1 : -N+1]


# In[29]:


from matplotlib.pyplot import legend

t = np.arange(N-1, len(bhp_returns))
plot(t, bhp_returns[N-1:], lw=1.0, label='bhp returns')
plot(t, smooth_bhp, lw=2.0, label='smooth bhp')
plot(t, vale_returns[N-1:], lw=1.0, label='vale returns')
plot(t, smooth_vale, lw=2.0, label='smooth vale')
legend(loc='best')
show()


# 图中折线有交叉，这些**交叉点可能是股价趋势的转折点**，至少可以表明BHP和VALE之间的股价关系发生了变化。这些转折点可能会经常出现，我们可以利用他们预测未来的股价走势。

# In[30]:


import matplotlib.pyplot as plt

# 使用多项式拟合平滑后的数据
K = 5
t = np.arange(N-1, len(bhp_returns))
poly_bhp = np.polyfit(t, smooth_bhp, K)
poly_vale = np.polyfit(t, smooth_vale, K)

fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.plot(t, smooth_bhp, label="smooth bhp")
poly_bhp_value = np.polyval(poly_bhp, t)
ax1.plot(t, poly_bhp_value, label='poly bhp')
plt.legend()

ax2 = fig.add_subplot(212)
ax2.plot(t, smooth_vale, label="smooth vale")
poly_vale_value = np.polyval(poly_vale, t)
ax2.plot(t, poly_vale_value, label='poly vale')
plt.legend()
show()


# In[31]:


# 得到交叉点的x坐标
# 通过求多项式函数差，再求根
poly_sub = np.polysub(poly_bhp, poly_vale)
xpoints = np.roots(poly_sub)
print "Intersection points:", xpoints


# In[32]:


# 判断是否为实数
# select选出实数
reals = np.isreal(xpoints)
print "Real number?",reals

xpoints = np.select([reals], [xpoints])
xpoints = xpoints.real
print "Real intersection points:", xpoints


# In[33]:


# 去除0元素
# trim_zeros函数可以去掉一维数组中开头和末尾为0的元素
print "Sans 0s", np.trim_zeros(xpoints)


# In[ ]: